{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "1ZlxP6qNoohn"
   },
   "source": [
    "# Text-Based Ideal Points using NumPyro\n",
    "\n",
    "## ___Szymon Sacher & Keyon Vafa___ This notebook replicates Text-Based Ideal Point model (Vafa, Naidu & Blei, 2020)\n",
    "\n",
    "This notebook is designed to run on Google Colab.\n",
    "\n",
    "**IMPORTANT:** To save this code and your results, make sure you copy to your personal Google Drive. Under \"File\", select \"Save a copy in Drive\".\n",
    "\n",
    "Use this Colab notebook to run a NumPyro implementation of the [text-based ideal point model (TBIP)](https://www.aclweb.org/anthology/2020.acl-main.475/) on a corpus of political text. The [Github repository is more complete](https://github.com/keyonvafa/tbip).\n",
    "\n",
    "See also [Tensorflow implementation](https://colab.research.google.com/drive/1_KkVI2lGtPdgsHSKDIMhSLCKkHvBQ4LO#scrollTo=BLgkfNNZUZv5) on which this notebook is based."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "mV0fFe6HpZXb"
   },
   "source": [
    "The [TBIP](https://www.aclweb.org/anthology/2020.acl-main.475/) is an unsupervised probabilistic topic model that analyzes texts to quantify the political positions of its authors. The model does not use political parties or votes, nor does it require any text labelled by ideology. Given a corpus of political text and the authors of each document, the TBIP estimates the latent political positions of the authors of texts and how per-topic word choice changes as a function of the political position of the author (\"ideological topics\"). [Refer to the paper for more information](https://www.aclweb.org/anthology/2020.acl-main.475/).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "bkhCDTQrprUf"
   },
   "source": [
    "## Getting started\n",
    "\n",
    "First, **make sure you are running this Colab using a GPU**. Go to the \"Runtime\" menu, and click \"Change runtime type\". If the \"Hardware accelerator\" is listed as \"None\" or \"TPU\", change to \"GPU\". Click \"Save\" and you're ready to go. Also, as described in the first cell, make sure this code is copied to your personal Google Drive."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "JPs53S8AqFzu"
   },
   "source": [
    "## Install NumPyro \n",
    "\n",
    "[NumPyro](https://github.com/pyro-ppl/numpyro) is a probabilistic programming framework powered by JAX for autograd and JIT compilation to GPU/TPU/CPU.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "iMkS4Z0bqFQx"
   },
   "outputs": [],
   "source": [
    "%%capture\n",
    "%pip install -q numpyro@git+https://github.com/pyro-ppl/numpyro\n",
    "%pip install optax"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "7WFN02jWqlhl"
   },
   "source": [
    "## Clone TBIP repository\n",
    "\n",
    "Below we clone the [Github repo for the TBIP](https://github.com/keyonvafa/tbip). This is where the data resides."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "NqVYM2QLqkhy",
    "outputId": "5ba84573-59f4-49fe-a9d0-84db135f31c4"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fatal: destination path 'tbip' already exists and is not an empty directory.\n"
     ]
    }
   ],
   "source": [
    "! git clone https://github.com/keyonvafa/tbip"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WAzgvHSbqxuJ"
   },
   "source": [
    "## Hyperparameters and Initialization\n",
    "\n",
    "We start setting some hyperparameters. We fix the number of topics $K = 50$. We also set a random seed for reproducibility."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "id": "TfFvLHY1q-yp"
   },
   "outputs": [],
   "source": [
    "from jax import random\n",
    "\n",
    "num_topics = 50\n",
    "rng_seed = random.PRNGKey(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "K8bzToM9rJxl"
   },
   "source": [
    "The next cell provides the data directory. The directory in the cell below links to speeches from the 114th Senate session from the `tbip` repo.\n",
    "\n",
    "To use your own corpus, upload the following four files to the Colab working directory:\n",
    "\n",
    "* `counts.npz`: a `[num_documents, num_words]` [sparse CSR matrix](https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.sparse.csr_matrix.html) containing the word counts for each document.\n",
    "* `author_indices.npy`: a `[num_documents]` vector where each entry is an integer in the set `{0, 1, ..., num_authors - 1}`, indicating the author of the corresponding document in `counts.npz`.\n",
    "* `vocabulary.txt`: a `[num_words]`-length file where each line denotes the corresponding word in the vocabulary.\n",
    "* `author_map.txt`: a `[num_authors]`-length file where each line denotes the name of an author in the corpus.\n",
    "\n",
    "See [Senate speech clean data](https://github.com/keyonvafa/tbip/tree/master/data/senate-speeches-114/clean) for an example of what the four files look like for Senate speeches. [Our setup script](https://github.com/keyonvafa/tbip/blob/master/setup/senate_speeches_to_bag_of_words.py) \n",
    "contains example code for creating the four files from unprocessed data for Senate speeches.\n",
    "\n",
    "**IMPORTANT:** If you are using your own corpus, change the following line to `data_dir = '.'` after uploading the four files to the Colab working directory.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "id": "X7M-44MgrONy"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import sparse\n",
    "\n",
    "import jax\n",
    "import jax.numpy as jnp\n",
    "\n",
    "dataPath = \"tbip/data/senate-speeches-114/clean/\"\n",
    "\n",
    "# Load data\n",
    "author_indices = jax.device_put(\n",
    "    jnp.load(dataPath + \"author_indices.npy\"), jax.devices(\"gpu\")[0]\n",
    ")\n",
    "\n",
    "counts = sparse.load_npz(dataPath + \"counts.npz\")\n",
    "\n",
    "with open(dataPath + \"vocabulary.txt\", \"r\") as f:\n",
    "    vocabulary = f.readlines()\n",
    "\n",
    "with open(dataPath + \"author_map.txt\", \"r\") as f:\n",
    "    author_map = f.readlines()\n",
    "\n",
    "author_map = np.array(author_map)\n",
    "\n",
    "num_authors = int(author_indices.max() + 1)\n",
    "num_documents, num_words = counts.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Si3avI3CrTBf"
   },
   "source": [
    "[In the paper](https://www.aclweb.org/anthology/2020.acl-main.475/), the parameters are pre-initialized with [Poisson factorization](https://arxiv.org/abs/1311.1704). Most of the time, we find this doesn't make a big difference for the learned ideal points, but it helps to interpret the ideological topics. \n",
    "\n",
    "Below, we initialize with Scikit-Learn's non-negative matrix factorization (NMF) implementation. Although we find that Poisson factorization learns more interpretable topics, we use Scikit-Learn's NMF implementation here because it is faster. To use Poisson factorization, see the [code in the Github repo](https://github.com/keyonvafa/tbip/blob/master/setup/poisson_factorization.py). \n",
    "\n",
    "If you would like to skip this pre-initialization step, set `pre_initialize_parameters = False` in the cell below. (Pre-initialization is recommended.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "id": "CEco12NLrkw4"
   },
   "outputs": [],
   "source": [
    "pre_initialize_parameters = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "iPpoenFXrtcJ"
   },
   "source": [
    "If you are pre-initializing parameters, the following cell might take a minute or so."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "id": "I-jdavQ7rpXM"
   },
   "outputs": [],
   "source": [
    "# Fit NMF to be used as initialization for TBIP\n",
    "from sklearn.decomposition import NMF\n",
    "\n",
    "if pre_initialize_parameters:\n",
    "    nmf_model = NMF(\n",
    "        n_components=num_topics, init=\"random\", random_state=0, max_iter=500\n",
    "    )\n",
    "    # Define initialization arrays\n",
    "    initial_document_loc = jnp.log(\n",
    "        jnp.array(np.float32(nmf_model.fit_transform(counts) + 1e-2))\n",
    "    )\n",
    "    initial_objective_topic_loc = jnp.log(\n",
    "        jnp.array(np.float32(nmf_model.components_ + 1e-2))\n",
    "    )\n",
    "else:\n",
    "    rng1, rng2 = random.split(rng_seed, 2)\n",
    "    initial_document_loc = random.normal(rng1, shape=(num_documents, num_topics))\n",
    "    initial_objective_topic_loc = random.normal(rng2, shape=(num_topics, num_words))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "id": "n8SUSFpLuXRP"
   },
   "source": [
    "## Perform Inference\n",
    "\n",
    "We perform inference using [variational inference](https://arxiv.org/abs/1601.00670) with [reparameterization](https://arxiv.org/abs/1312.6114) [gradients](https://arxiv.org/abs/1401.4082). We provide a brief summary below, but encourage readers to [refer to the original paper](https://www.aclweb.org/anthology/2020.acl-main.475/) for a more complete overview.\n",
    "\n",
    "It is intractable to evaluate the posterior distribution $p(\\theta, \\beta, \\eta, x | y)$, so we approximate the posterior with a distribution $q_\\phi(\\theta, \\beta,\\eta,x)$, parameterized by $\\phi$. How do we set the values $\\phi$? We want to minimize the KL-Divergence between $q$ and the posterior, which is equivalent to maximizing the ELBO:\n",
    "$$\\mathbb{E}_{q_\\phi}[\\log p(y, \\theta, \\beta, \\eta, x) - \\log q_{\\phi}(\\theta, \\beta, \\eta, x)].$$\n",
    "We set the variational family to be the mean-field family, meaning the latent variables factorize over documents $d$, topics $k$, and authors $s$:\n",
    "$$q_\\phi(\\theta, \\beta, \\eta, x) = \\prod_{d,k,s} q(\\theta_d)q(\\beta_k)q(\\eta_k)q(x_s).$$\n",
    "We use lognormal factors for the positive variables and Gaussian factors for the real variables:\n",
    "$$q(\\theta_{dk}) = \\text{LogNormal}(\\mu_{\\theta_{dk}}\\sigma^2_{\\theta_{dk}})$$\n",
    "$$q(\\beta{dv}) = \\text{LogNormal}(\\mu_{\\beta_{kv}}, \\sigma^2_{\\beta_{kv}})$$\n",
    "$$q(\\eta_{kv}) = \\text{Normal}(\\mu_{\\eta_{kv}}, \\sigma^2_{\\eta_{kv}})$$\n",
    "$$q(x_s) = \\text{Normal}(\\mu_{x_s}, \\sigma^2_{x_s}).$$\n",
    "\n",
    "Thus, our goal is to maximize the ELBO with respect to $\\phi = \\{\\mu_\\theta, \\sigma_\\theta, \\mu_\\beta, \\sigma_\\beta,\\mu_\\eta, \\sigma_\\eta, \\mu_x, \\sigma_x\\}$. \n",
    "\n",
    "In the cell below, we define the model and the variational family (guide)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "id": "iRayId2goohr"
   },
   "outputs": [],
   "source": [
    "from numpyro import param, plate, sample\n",
    "import numpyro.distributions as dist\n",
    "from numpyro.distributions import constraints\n",
    "\n",
    "# Define the model and variational family\n",
    "\n",
    "\n",
    "class TBIP:\n",
    "    def __init__(self, N, D, K, V, batch_size, init_mu_theta=None, init_mu_beta=None):\n",
    "        self.N = N  # number of people\n",
    "        self.D = D  # number of documents\n",
    "        self.K = K  # number of topics\n",
    "        self.V = V  # number of words in vocabulary\n",
    "        self.batch_size = batch_size  # number of documents in a batch\n",
    "\n",
    "        if init_mu_theta is None:\n",
    "            init_mu_theta = jnp.zeros([D, K])\n",
    "        else:\n",
    "            self.init_mu_theta = init_mu_theta\n",
    "\n",
    "        if init_mu_beta is None:\n",
    "            init_mu_beta = jnp.zeros([K, V])\n",
    "        else:\n",
    "            self.init_mu_beta = init_mu_beta\n",
    "\n",
    "    def model(self, Y_batch, d_batch, i_batch):\n",
    "        with plate(\"i\", self.N):\n",
    "            # Sample the per-unit latent variables (ideal points)\n",
    "            x = sample(\"x\", dist.Normal())\n",
    "\n",
    "        with plate(\"k\", size=self.K, dim=-2):\n",
    "            with plate(\"k_v\", size=self.V, dim=-1):\n",
    "                beta = sample(\"beta\", dist.Gamma(0.3, 0.3))\n",
    "                eta = sample(\"eta\", dist.Normal())\n",
    "\n",
    "        with plate(\"d\", size=self.D, subsample_size=self.batch_size, dim=-2):\n",
    "            with plate(\"d_k\", size=self.K, dim=-1):\n",
    "                # Sample document-level latent variables (topic intensities)\n",
    "                theta = sample(\"theta\", dist.Gamma(0.3, 0.3))\n",
    "\n",
    "            # Compute Poisson rates for each word\n",
    "            P = jnp.sum(\n",
    "                jnp.expand_dims(theta, 2)\n",
    "                * jnp.expand_dims(beta, 0)\n",
    "                * jnp.exp(\n",
    "                    jnp.expand_dims(x[i_batch], (1, 2)) * jnp.expand_dims(eta, 0)\n",
    "                ),\n",
    "                1,\n",
    "            )\n",
    "\n",
    "            with plate(\"v\", size=self.V, dim=-1):\n",
    "                # Sample observed words\n",
    "                sample(\"Y_batch\", dist.Poisson(P), obs=Y_batch)\n",
    "\n",
    "    def guide(self, Y_batch, d_batch, i_batch):\n",
    "        # This defines variational family. Notice that each of the latent variables\n",
    "        # defined in the sample statements in the model above has a corresponding\n",
    "        # sample statement in the guide. The guide is responsible for providing\n",
    "        # variational parameters for each of these latent variables.\n",
    "\n",
    "        # Also notice it is required that model and the guide have the same call.\n",
    "\n",
    "        mu_x = param(\n",
    "            \"mu_x\", init_value=-1 + 2 * random.uniform(random.PRNGKey(1), (self.N,))\n",
    "        )\n",
    "        sigma_x = param(\n",
    "            \"sigma_y\", init_value=jnp.ones([self.N]), constraint=constraints.positive\n",
    "        )\n",
    "\n",
    "        mu_eta = param(\n",
    "            \"mu_eta\", init_value=random.normal(random.PRNGKey(2), (self.K, self.V))\n",
    "        )\n",
    "        sigma_eta = param(\n",
    "            \"sigma_eta\",\n",
    "            init_value=jnp.ones([self.K, self.V]),\n",
    "            constraint=constraints.positive,\n",
    "        )\n",
    "\n",
    "        mu_theta = param(\"mu_theta\", init_value=self.init_mu_theta)\n",
    "        sigma_theta = param(\n",
    "            \"sigma_theta\",\n",
    "            init_value=jnp.ones([self.D, self.K]),\n",
    "            constraint=constraints.positive,\n",
    "        )\n",
    "\n",
    "        mu_beta = param(\"mu_beta\", init_value=self.init_mu_beta)\n",
    "        sigma_beta = param(\n",
    "            \"sigma_beta\",\n",
    "            init_value=jnp.ones([self.K, self.V]),\n",
    "            constraint=constraints.positive,\n",
    "        )\n",
    "\n",
    "        with plate(\"i\", self.N):\n",
    "            sample(\"x\", dist.Normal(mu_x, sigma_x))\n",
    "\n",
    "        with plate(\"k\", size=self.K, dim=-2):\n",
    "            with plate(\"k_v\", size=self.V, dim=-1):\n",
    "                sample(\"beta\", dist.LogNormal(mu_beta, sigma_beta))\n",
    "                sample(\"eta\", dist.Normal(mu_eta, sigma_eta))\n",
    "\n",
    "        with plate(\"d\", size=self.D, subsample_size=self.batch_size, dim=-2):\n",
    "            with plate(\"d_k\", size=self.K, dim=-1):\n",
    "                sample(\"theta\", dist.LogNormal(mu_theta[d_batch], sigma_theta[d_batch]))\n",
    "\n",
    "    def get_batch(self, rng, Y, author_indices):\n",
    "        # Helper functions to obtain a batch of data, convert from scipy.sparse\n",
    "        # to jax.numpy.array and move to gpu\n",
    "\n",
    "        D_batch = random.choice(rng, jnp.arange(self.D), shape=(self.batch_size,))\n",
    "        Y_batch = jax.device_put(jnp.array(Y[D_batch].toarray()), jax.devices(\"gpu\")[0])\n",
    "        D_batch = jax.device_put(D_batch, jax.devices(\"gpu\")[0])\n",
    "        I_batch = author_indices[D_batch]\n",
    "        return Y_batch, I_batch, D_batch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "5ImzSMzSviB3"
   },
   "source": [
    "## Initialization\n",
    "\n",
    "Below we initialize an instance of the TBIP model, and associated SVI object. The latter is used to compute the Evidence Lower Bound (ELBO) given current value of guide's parameters and current batch of data.\n",
    "\n",
    "We optimize the model using Adam with exponential decay of learning rate.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "id": "LxGTUSTXoohs"
   },
   "outputs": [],
   "source": [
    "# Initialize the model\n",
    "from jax import jit\n",
    "from optax import adam, exponential_decay\n",
    "\n",
    "from numpyro.infer import SVI, TraceMeanField_ELBO\n",
    "\n",
    "num_steps = 50000\n",
    "batch_size = 512  # Large batches are recommended\n",
    "learning_rate = 0.01\n",
    "decay_rate = 0.01\n",
    "\n",
    "tbip = TBIP(\n",
    "    N=num_authors,\n",
    "    D=num_documents,\n",
    "    K=num_topics,\n",
    "    V=num_words,\n",
    "    batch_size=batch_size,\n",
    "    init_mu_theta=initial_document_loc,\n",
    "    init_mu_beta=initial_objective_topic_loc,\n",
    ")\n",
    "\n",
    "svi_batch = SVI(\n",
    "    model=tbip.model,\n",
    "    guide=tbip.guide,\n",
    "    optim=adam(exponential_decay(learning_rate, num_steps, decay_rate)),\n",
    "    loss=TraceMeanField_ELBO(),\n",
    ")\n",
    "\n",
    "# Compile update function for faster training\n",
    "svi_batch_update = jit(svi_batch.update)\n",
    "\n",
    "# Get initial batch. This informs the dimension of arrays and ensures they are\n",
    "# consistent with dimensions (N, D, K, V) defined above.\n",
    "Y_batch, I_batch, D_batch = tbip.get_batch(random.PRNGKey(1), counts, author_indices)\n",
    "\n",
    "# Initialize the parameters using initial batch\n",
    "svi_state = svi_batch.init(\n",
    "    random.PRNGKey(0), Y_batch=Y_batch, d_batch=D_batch, i_batch=I_batch\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "id": "RQfTWq9sxmxw"
   },
   "outputs": [],
   "source": [
    "# @title Run this cell to create helper function for printing topics\n",
    "\n",
    "\n",
    "def get_topics(\n",
    "    neutral_mean, negative_mean, positive_mean, vocabulary, print_to_terminal=True\n",
    "):\n",
    "    num_topics, num_words = neutral_mean.shape\n",
    "    words_per_topic = 10\n",
    "    top_neutral_words = np.argsort(-neutral_mean, axis=1)\n",
    "    top_negative_words = np.argsort(-negative_mean, axis=1)\n",
    "    top_positive_words = np.argsort(-positive_mean, axis=1)\n",
    "    topic_strings = []\n",
    "    for topic_idx in range(num_topics):\n",
    "        neutral_start_string = \"Neutral  {}:\".format(topic_idx)\n",
    "        neutral_row = [\n",
    "            vocabulary[word] for word in top_neutral_words[topic_idx, :words_per_topic]\n",
    "        ]\n",
    "        neutral_row_string = \", \".join(neutral_row)\n",
    "        neutral_string = \" \".join([neutral_start_string, neutral_row_string])\n",
    "\n",
    "        positive_start_string = \"Positive {}:\".format(topic_idx)\n",
    "        positive_row = [\n",
    "            vocabulary[word] for word in top_positive_words[topic_idx, :words_per_topic]\n",
    "        ]\n",
    "        positive_row_string = \", \".join(positive_row)\n",
    "        positive_string = \" \".join([positive_start_string, positive_row_string])\n",
    "\n",
    "        negative_start_string = \"Negative {}:\".format(topic_idx)\n",
    "        negative_row = [\n",
    "            vocabulary[word] for word in top_negative_words[topic_idx, :words_per_topic]\n",
    "        ]\n",
    "        negative_row_string = \", \".join(negative_row)\n",
    "        negative_string = \" \".join([negative_start_string, negative_row_string])\n",
    "\n",
    "        if print_to_terminal:\n",
    "            topic_strings.append(negative_string)\n",
    "            topic_strings.append(neutral_string)\n",
    "            topic_strings.append(positive_string)\n",
    "            topic_strings.append(\"==========\")\n",
    "        else:\n",
    "            topic_strings.append(\n",
    "                \"  \\n\".join([negative_string, neutral_string, positive_string])\n",
    "            )\n",
    "\n",
    "    if print_to_terminal:\n",
    "        all_topics = \"{}\\n\".format(np.array(topic_strings))\n",
    "    else:\n",
    "        all_topics = np.array(topic_strings)\n",
    "    return all_topics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "vBzuncDP0NN_"
   },
   "source": [
    "## Execute Training\n",
    "The code above was creating the model; below we actually run training. You can adjust the number of steps to train (`num_steps`, defined above) and the frequency at which to print the ELBO (`print_steps`, defined below).\n",
    "\n",
    "Here, we run our training loop. Topic summaries and ordered ideal points will print every 2500 steps. Typically in our experiments it takes 15,000 steps or so to begin seeing sensible results, but of course this depends on the corpus. These sensible results should be reached within a half hour. For the default corpus of Senate speeches, it should take less than 2 hours to complete 30,000 training steps (which is the default `num_steps`). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "xoTt8vn8ooht",
    "outputId": "eec708bc-cbd5-4d06-945c-79b06e4b6dc2"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Init loss: 14323.9385; Avg loss (last 100 iter):   953.5815:   5%|▍         | 2499/50000 [03:47<1:11:33, 11.06it/s]/usr/local/lib/python3.7/dist-packages/jax/_src/numpy/lax_numpy.py:3327: UserWarning: 'kind' argument to argsort is ignored; only 'stable' sorts are supported.\n",
      "  warnings.warn(\"'kind' argument to argsort is ignored; only 'stable' sorts \"\n",
      "Init loss: 14323.9385; Avg loss (last 100 iter):   634.0942: 100%|██████████| 50000/50000 [1:17:57<00:00, 10.69it/s]\n"
     ]
    }
   ],
   "source": [
    "# Run SVI\n",
    "import pandas as pd\n",
    "from tqdm import tqdm\n",
    "\n",
    "print_steps = 100\n",
    "print_intermediate_results = False\n",
    "\n",
    "rngs = random.split(random.PRNGKey(2), num_steps)\n",
    "losses = []\n",
    "pbar = tqdm(range(num_steps))\n",
    "\n",
    "\n",
    "for step in pbar:\n",
    "    Y_batch, I_batch, D_batch = tbip.get_batch(rngs[step], counts, author_indices)\n",
    "    svi_state, loss = svi_batch_update(\n",
    "        svi_state, Y_batch=Y_batch, d_batch=D_batch, i_batch=I_batch\n",
    "    )\n",
    "\n",
    "    loss = loss / counts.shape[0]\n",
    "    losses.append(loss)\n",
    "    if step % print_steps == 0 or step == num_steps - 1:\n",
    "        pbar.set_description(\n",
    "            \"Init loss: \"\n",
    "            + \"{:10.4f}\".format(jnp.array(losses[0]))\n",
    "            + f\"; Avg loss (last {print_steps} iter): \"\n",
    "            + \"{:10.4f}\".format(jnp.array(losses[-100:]).mean())\n",
    "        )\n",
    "\n",
    "    if (step + 1) % 2500 == 0 or step == num_steps - 1:\n",
    "        # Save intermediate results\n",
    "        estimated_params = svi_batch.get_params(svi_state)\n",
    "\n",
    "        neutral_mean = (\n",
    "            estimated_params[\"mu_beta\"] + estimated_params[\"sigma_beta\"] ** 2 / 2\n",
    "        )\n",
    "\n",
    "        positive_mean = (\n",
    "            estimated_params[\"mu_beta\"]\n",
    "            + estimated_params[\"mu_eta\"]\n",
    "            + (estimated_params[\"sigma_beta\"] ** 2 + estimated_params[\"sigma_eta\"] ** 2)\n",
    "            / 2\n",
    "        )\n",
    "\n",
    "        negative_mean = (\n",
    "            estimated_params[\"mu_beta\"]\n",
    "            - estimated_params[\"mu_eta\"]\n",
    "            + (estimated_params[\"sigma_beta\"] ** 2 + estimated_params[\"sigma_eta\"] ** 2)\n",
    "            / 2\n",
    "        )\n",
    "\n",
    "        np.save(\"neutral_topic_mean.npy\", neutral_mean)\n",
    "        np.save(\"negative_topic_mean.npy\", positive_mean)\n",
    "        np.save(\"positive_topic_mean.npy\", negative_mean)\n",
    "\n",
    "        topics = get_topics(neutral_mean, positive_mean, negative_mean, vocabulary)\n",
    "\n",
    "        with open(\"topics.txt\", \"w\") as f:\n",
    "            print(topics, file=f)\n",
    "\n",
    "        authors = pd.DataFrame(\n",
    "            {\"name\": author_map, \"ideal_point\": np.array(estimated_params[\"mu_x\"])}\n",
    "        )\n",
    "        authors.to_csv(\"authors.csv\")\n",
    "\n",
    "        if print_intermediate_results:\n",
    "            print(f\"Results after {step} steps.\")\n",
    "            print(topics)\n",
    "            sorted_authors = \"Authors sorted by their ideal points: \" + \",\".join(\n",
    "                list(authors.sort_values(\"ideal_point\")[\"name\"])\n",
    "            )\n",
    "            print(sorted_authors.replace(\"\\n\", \" \"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "id": "XL9A29cFooht"
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "neutral_topic_mean = np.load(\"neutral_topic_mean.npy\")\n",
    "negative_topic_mean = np.load(\"negative_topic_mean.npy\")\n",
    "positive_topic_mean = np.load(\"positive_topic_mean.npy\")\n",
    "authors = pd.read_csv(\"authors.csv\")\n",
    "authors[\"name\"] = authors[\"name\"].str.replace(\"\\n\", \"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "o-Om5Qk31ENx"
   },
   "source": [
    "For example, here is a graph of the learned ideal points. We don't label each point because there are too many to plot. Below we select some authors to label."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 183
    },
    "id": "1J8Sw4-7oohu",
    "outputId": "dbf8acf0-db94-4814-c6be-c370b6b116f5"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA8wAAACmCAYAAAAPrwbPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd1xUV9rA8d8wM/TO0HuR3osggtgQgkajJllLNKtmE9e8aWuSTdcku9lN7GXtUWNvUWMvqKBGrNiwoSiIiCCCSm8z7x9+5i6oyTYjxfP9Rx2H4d47Z+6c55znPEem0Wg0CIIgCIIgCIIgCILQjE5LH4AgCIIgCIIgCIIgtEYiYBYEQRAEQRAEQRCExxABsyAIgiAIgiAIgiA8hgiYBUEQBEEQBEEQBOExRMAsCIIgCIIgCIIgCI8hAmZBEARBEARBEARBeAwRMAuCIAiCIAiCIAjCY4iAWRAEQRAEQRAEQRAeQwTMgiAIgiAIrcSZM2eYPHky165dA0Cj0bTwEQlC25Gdnc25c+cA8dkRnhwRMAuCIAiCILQSZmZmXL9+nbfffpudO3cik8la+pAEoU0oKipi2bJljBkzhvT0dCorKwFQq9UtfGRCWyfTiOEXQRAEQRCEp06j0SCTyVCr1ejo6NDY2IhcLgdg4cKFrFy5kpdffpnXXntNBM6C8G+aN28eGRkZODo68pe//KWlD0doB+Tjx48f39IHIQiCIAiC8KzQzlVog2Dtnzo6OlIQHRYWRmVlJbt376aiooLw8PAWO15BaE0e/vxoNTY2oqOjQ0REBMbGxixevJj79+8TExPTEocptCOKlj4AQRAEQRCEZ4V2NhngxIkTpKenU19fT319Pc8//zx+fn7o6uoCMGjQIGpqali4cCHJyck4Ojq25KELQotr+vm5dOkSt27dQqlU4ubmhq2trfS8nj17Ul1dzZ///GdCQkLo2rWryNJopyoqKjA2Nm7WNp40kZItCIIgCILwFN24cYPPP/+cy5cvk5CQwK1btygqKqKgoICRI0fy1ltvSc+9fPkyX3zxBY6OjkycOLEFj1oQWo428wKgoKCAr7/+mvz8fAwNDTl//jwWFhZ0796dr776qtnPjR49mjt37jB37lwsLS1b4tCF34harWbDhg1s3ryZcePG4erq2ixL50kSRb8EQRAEQRCekqlTp9KzZ0/s7e1ZuXIl48aN4/vvv2fz5s0kJCSwevVqVq5cKT3f3d2dXr16ce3aNXJyclrwyAXh6Xs4/fqrr74iKSkJIyMjJk2axHfffcehQ4dITk5m7dq1zJw5s1mRr08//ZSsrCwuXbrU7PWEtqnp+6ejo0OnTp1wcnLiyy+/ZP78+cCjqfpPgljDLAiCIAiC8Bs5efIkNTU1WFhYkJWVxddff02vXr349ttvMTMzQ0dHB5lMhkwmw9/fn6ysLDIyMujduzf6+vro6OhQU1PDkSNH8PPzw9XVtaVPSRCemrKyMgwMDKioqGDKlCksXbqUjRs3MmTIEFQqFRYWFujr6xMcHIxGo2Hp0qV07twZW1tb1Go15ubmZGVlce7cOVJSUkRadhunff9ycnKwtLTExMSE+Ph4AGbMmIGenh5BQUFihlkQBEEQBKEtyM3N5e9//zv6+voABAYGMmDAAIqKikhPT3/k+S4uLnTr1o2ysjLS0tKkxzt16sS9e/e4ffs2IGbJhGfDihUr+OSTTwAwNjYmOjoaKysrTp8+DUB9fb30XHNzc373u99hbm7O6tWrpcc1Gg3dunWjrq6O0tLSp3sCwhNz584dAOrq6li8eDGvvvoqRUVFACiVSgYOHMj777/P5MmT2bVr1xPfSkwEzIIgCIIgCL8BFxcX7ty5w5UrV6THXnnlFerq6khNTeXu3bvIZDI0Go3UwUtMTKSiooKKigrgQQcRICEhgfPnzwO/TcqhILQ2RUVFBAUFSf+OiIjg+eefZ8qUKdTW1qJUKpsFRtbW1sTFxXHw4EEqKyul7A1zc3NKSkqwsLBoidMQ/kfHjx9n7Nix1NbWoquri7u7O46OjlIKtnYA8eWXXyYmJobly5dz+fLlJ3oMImAWBEEQBEH4Ddy7dw9fX19KSkqABx07Z2dnkpKSOHfuHPv27QMeBMDa6q56enoYGhpKr6GtmN3Q0ECHDh2e8hkIwtOnDYJNTU3ZtWuX9LiJiQl9+vTBxMSEb7/99pGfMzAwQE9PDzMzM5RKJY2NjQCEhoZSVVXFjRs3ns4JCE/UxYsXMTMzQ09PD3iQcZOQkMD+/fs5ffo0MplMGlj8+OOPyc7O5sSJE0/0GETALAiCIAiC8BuwsLDg3r17nDp1CngQ9MKD7aIsLCzYt28feXl5wIM9ZNVqNatXr8bY2Jhu3bpJjwP06dOH4ODgFjgLQXi6tINHvr6+6Onpce7cOen/vL29GTx4MOvXrycnJwcdHR3pc1VaWsqxY8fw9PREoVAgl8uBB6nbr776KjY2Nk//ZIT/mnbm2NfXl9OnT1NZWQk8GETs1q0bLi4uzJkzR3qsoaEBV1dXunXrxqZNm5q9xv9KBMyCIAjCM0Ws/xSeBm2gO2DAADZt2iSlkNbX12NoaMiLL75IXl4eO3fuBEAul7Nt2zZSU1MZPnw4Dg4OaDQaqdMfHh6Ov79/i52PIDxthoaG1NfXk52dLT2mDZYCAwP561//CoBCoQBgy5YtGBoa8uabbzbbj9fJyYmBAwdKM5RC26BdeiKTyXB2dpYGHgH8/Pzo2bMnV69eZevWrcA/MxOGDBlCTk4Oubm5T2z5igiYBUEQhGdGbW2tWP8pPFHawFhL22nTBrodO3bEwsKCyZMnA//sBD733HN4e3tz4sQJtm7dyujRoxk3bhz9+vVj2LBhzZ4rCM+isLAwTE1N2b17t1TwDh7UBhg6dCgnT54kIyOD8vJyRo0axYIFCxgyZAienp6PvJZ2aYPQ9ri5uXHr1i2uXLnSbM16fHw8wcHB/PDDD9TV1UnvsYmJCSEhIeTn5z+xYxABsyAIgvBMqKysZPjw4fz1r38lKysL4IlX0hSeHQ8Hxj/88ANVVVXNZrbgQSGiYcOGsXz5cnJyclAoFNLPDhs2jPPnzzN27Fj09PRIT0/nlVdeAUQmhNC+PTzQ9Ev//9Zbb7F//34OHjwoVcWWyWR07NiRxMRERowYQUxMDMbGxmzbto0+ffr85scuPD1qtRorKyvi4uJYtWoVZWVl0v85OjqSmJhIdXU1ixYtkh53cXEhPz9fujc/CWIfZkEQBOGZoKurS0BAAGfOnGHRokUEBgZib2/f0ocltFHa2d/bt28zfvx41q9fL+3/2pRcLsfZ2ZkrV66wZs0a4uPjMTc3B8DW1hZLS0vee+89hgwZgq6uLo2NjVJ1X0FojxobG6VgpqGhQRpkUqvVUrvX0dFBo9Hg6OhIbm4uBw8exMrKSpo9NjQ0RFdXF319fb7++muGDx+Orq5us9cQWj+1Wo1Go/nV90wmkxEQEMCsWbMwNTUlODhYajPW1tbcvn2bdevW0bNnT2lv+6KiIqysrPDw8HgixykCZkEQBKHdeVynSaPRYGNjQ8eOHbl58yZz5szBwcHhsel7gvCvVFVVMWvWLH788UcMDAyYO3cuLi4uj32uoaEh3bt3Z926dVy+fBlzc3OcnJyABwVtLC0tpVnnh2eoBaG90AbHOjo6XL16lXHjxpGWlkZZWRmBgYHSFmvae7f271FRUZw8eZIDBw7g5uaGg4MDAA4ODiQmJmJtbS0+P22QduBEJpNRUFBAVlYWzs7OzdqAtk0YGRmhUCiYO3cuHTp0kAJh7br00tJSQkNDUalUwIN9uSMiIp5YexABsyAIgtBuaANlmUxGY2MjR44cQa1WY2hoKKXC6unp0aVLF06cOEFmZiY2NjY4Ozu39KELrZh21vfhx3788UcOHTpEREQEPXr0aDZb9vBz9fT0CA8Pp7CwkOnTp+Pt7Y2dnZ0006Ztt+2JWq3m6NGjWFlZoVQqW/pwhBamo6ODWq2muLiYV199FWdnZ+7fv8/atWuRyWREREQAzYs9qdVqDAwMpO3ZZsyYgaWlJV5eXlKb0mg0IiujDdJWOB8/fjx///vfKSoqIjY2ttm2evDP9hAeHk5WVhaZmZnIZDJ8fHyAB0XdUlJSpGAZHsw8P8nBExEwC4IgCO2CttMEsHjxYt5++21OnjzJ999/T3l5OSEhIejp6VFXV4dcLsfV1ZWff/6ZvLw8unbt+kTXO7Vla9euZdmyZYSGhj7ScXkWNW1XmZmZyOVyNBoNhoaGqFQqjh49SmNjI88///wvdtC0j6tUKmJjY9HX16eurg47O7t2fY23bNnCDz/8wO7duwkNDcXMzKylD0loQfn5+QwYMICcnBx69erFhx9+SNeuXTExMWHKlCn06dNHWq6gpQ2WzM3NiYuLQ6FQcO7cOWpra6V9yUWg3DZoswCaLmcZPXo0JSUlTJgwgV69emFpafnY72LtoGVoaCiVlZV89913ODk5YWVlhZGREcC/TO3+X4iAWfhFGo2G3NxcKW3R0tKypQ9JEAThF8lkMjIzMxk1ahTHjh3jo48+4rXXXqOqqoqtW7dibGxMcHCw9GVsbW1NaWkpJ0+exN7eHldX1xY+g9ZBqVSyb98+duzYgYmJCe7u7i19SC1KJpNx7NgxXn/9dXbv3s3WrVvZvXs33bp1w9PTk9LSUs6dO4etrS3u7u6/uoZS26ELDAzEz8+v3QbL2vP08fGhd+/ebN26lX379knt6bfs2AotT6PRoFarHxlAun//PlevXmX37t384Q9/wNHRET09PVxcXDh8+DAnTpz4xaJd2s9VSEgIcXFxeHl5ifTrNkJbwFCbBaCtbF5eXs6PP/7IunXrcHBwQKFQYGBgQEVFxSNVzbXvtYmJCREREZiZmXHy5ElKS0sJCQkBftuBE9HShF8kk8mwsLCguLiYDz74gNWrV7f0IQmCIDRz7949qeI1wMSJE7l37x4LFiwgJSUFe3t7hg0bRllZGXfv3pWepx3pTklJoaKighs3bjz1Y28NHq4S3tjYiLe3N19//TUdO3Zk7NixpKamStVpn0WXL19m/Pjx9OzZkzVr1rB06VIKCwv58MMPqaiooHfv3qhUKjZs2EBtba1UrOhxmqaaticPVzzWnl99fT26urp88cUXeHt78+c//5mcnJx2d/7CP2kDW7lcTnFxMampqVy8eBEAZ2dnXnrpJZRKJZcuXQIetB1LS0vefPNN9u/fz6FDh4BHq8Q3DY6VSqW097LQ+mmXm9TU1DBv3jzGjBnD2bNnqa+vp6CggG3btjFz5kwmTpzIwIEDGTRoED///POv7mIxePBgPv/8c1566aWncg5ihlmQPDzi29jYiKGhIcnJyejp6TFjxgyqqqoIDAwUa5EEQWgVPvroI6qqqoiMjAQerGVKS0tDpVJJo84zZ87k9OnTJCUlERwcLKXYNjY2YmpqyqlTp8jKyqJfv37PzMxX00rMFRUVVFVVoa+vL3VKjYyMiIyMJC8vj/T0dAwNDfH29m7X1+dx65QBlixZglqt5uuvv0ahULBs2TJSU1OJj48nJiYGW1tb7t+/T2ZmJgCBgYFA+wuKH0cbHGmv25IlS8jMzOT27dt4enpK2RxmZmZ07NiR/fv3c+DAAeLi4qQ0SqF90bb7SZMm8dFHH5GXl8fcuXMpLi7GxcWFoKAgiouLWblyJcOGDZP6k7a2tuTl5bF48WJGjRr1THx+niXr16/np59+Ii8vjy+//JIePXqgr6+Pnp4es2fPprGxERcXF+Lj4ykoKCAjI4OUlJRfjTcUCsVTi0dEwCyg0WiardEqLy9HT09P+rdGo8HPzw8dHR327NlDQUEBsbGxLXnIgiA847Qd9QsXLnDjxg0SExOBBwHz2bNnOXPmDNevX+eDDz7gxo0b2NjYAFBSUoK1tTXGxsZSyqCxsTH79u0jKSkJAwODljyt35w24NXe3ydMmMC0adPYunUrBQUFuLm5YWJiIq3zDgkJISMjgwsXLhAWFoaZmVm7DJqbfgdWVFSgVCqlc9y7dy+6urrIZDJeffVVioqK+Pvf/87QoUOltEFHR0cOHTpEZmYmPXv2bPftSEt7jdasWcPIkSO5efMmN27c4Pvvv8fIyAhfX1+USiUajQaFQkFUVBTfffcdrq6u0sCC0LZp+5BN7wkrV65ky5YtTJo0ibfffpuEhAQmTpxIaWkpXbt2xcbGhrS0NG7cuEGXLl1Qq9Xo6uri6OjI3bt3iYuLExMzbdTD65S11q5dy/r16zExMeGPf/wjGo0GAwMDoqOjGTRoEC+99BKhoaEEBQVx7949zp8/z8CBA1tNJoEImJ9x2g6jTCbj1KlTjB8/nq1bt3Lq1Cns7OxQqVTU19cjl8vp0KEDFRUVrFixgujo6Ef2mhQEQXhatF/Ghw8fpr6+ns6dO0tBT0BAAAsXLuTo0aOMGTOGyZMnM2DAANRqNWvWrGHFihUEBwdjZmaGrq4uWVlZXLhwgZdffrndr4nTXreNGzfyyiuvUF1dzciRI6mvr2fbtm3cv3+fLl26IJfLUavVGBsbo1AoOHLkiFRwpb0Fy/Dgupw9e5a33npLSiGNi4sDHhT62rBhA3v27GHs2LF8/PHHuLm5UVdXx5IlSzA1NcXR0RELCwt69+79TFVcLysrY/jw4Wzbto3PPvuMzz//nOTkZO7evcumTZtITk7GzMxMqlpvYWHB7du32bVrF4MHD26Xgy/PkqaZKocPHyYjI4OAgACmTp1Knz59SE5OJjMzkwkTJlBRUcHw4cPx9vbG1NQUmUzGkiVLSExMlGrk2NrakpiYKILlNqrpNlH5+fncv38fhUKBrq4u/v7+nDhxgvv379OrVy8MDQ2l9qOrq8v9+/fRaDTs2bOH+fPn88ILLxAVFdVq7g/tu2cg/Es6Ojrcv3+fMWPGMGzYMFxcXPDz8+PQoUN88sknANJG8EZGRiQnJxMUFMS0adNa+MgFQXiWaUexvb292bZtG3K5HIVCQUNDA46OjgwZMgQbGxu8vb0B0NfXJyUlhTlz5hAQEMBrr73GhAkTAIiJicHGxoaampoWO5/fWmlpqXTNli1bxkcffcTgwYNZs2YNvXv3Zvz48djY2HD69GlKS0uBfwbXPXv2xNLSkkuXLlFZWdli5/AkPbzmtqysjG+++YaQkBBCQkJYtmwZ3377LQDx8fGoVCoSExMZOHCgtO/nmTNn2LlzJ+fOnQOgc+fO+Pr6Pt0TeYo0Gg0NDQ188803UjuorKykvLyc+Ph4+vbti0wmw8jICAMDA27dutWsboA2PXvIkCFcv36dkydPSnusCm2L9j2Ty+XU1dVx9OhRXn/9dS5fvgxAXl4e+vr6fPnll7z22mv4+fmxceNGkpKSqK+vR19fnx49eqBSqZgyZcovvr7QelVUVDxS20Iul3P//n3efvttRowYwZgxYxg9ejRnzpzB2tqaF154AY1Gw86dO6XnA9y4cYPVq1czatQovvzySwYPHtzq0vLFDPMzSrtX5OXLlxk5ciQ3btwgNTWVxMREOnfuTHl5OVu2bCEkJESqHCuTyTAzM6OmpoYjR47g5+eHnZ1dC5+JIAjPIu0XqUqlYu3atejr6xMcHExDQwNyuZzQ0FDWrVtHWVkZQUFBGBsb09jYiIGBAd26dSM4OJhBgwZJWwTFxsa22y1vtCP2PXr0QC6XY2trS2pqKr6+vgQGBqKrq8uNGzdYuXIl+vr6DBs2TJo1amhoQKFQUFtby6pVqxg1alSb3n7r4TW3a9euJT8/n8uXL2Ntbc2f/vQnYmJicHNzY+LEiYSFhREdHU1ZWRmpqamkpaVRUlLC2rVrmTJlCr1795auV3un3RN3+vTp+Pr64uDggImJCdbW1syYMUPaEmjmzJksXLgQlUrFG2+88Uh6emNjI7m5uSiVSkJCQlpVp1j4ddqMAO17dvz4cd555x2uXr3KmDFjpGyVq1evMnPmTFQqFVOnTqVfv34YGBhQUFDAhg0bsLa2xsnJicjISOk+3JRoE61beXk5zz33HPfu3Wu2RLOyspJ3332Xuro6vvvuO5KTk7ly5QozZ87k+eefJzIykoMHD5KXl4e3tzdWVlbAgyJutbW12NvbM3HiRGk/7taUgdL+7/BCM+Xl5Xz99dfs378fADs7OwICAujQoUOzkeCioiL09fWlNBntFyWAj48Penp6FBQUPP0TEAThmfKvZhp0dXVJTk5m9erVVFdXo6urS11dHQqFgj/84Q8cOXKEw4cPA0jBsZ6eHp07d0ZXV1faT7c9bpunvXZmZmacP39eWgtma2vLyy+/zK5du8jMzOQf//gHvXv35ubNm3To0IE9e/ZQXV0N/LPjmpSUhL6+PkePHm2Zk3lCtIHtlStXeOGFF5g0aRIzZ86UMqq0UlJSCAkJYfbs2VRXV/P666/z9ddfo6enx4kTJ6ioqGDdunW89957Urt6FhQVFVFXVyd1dGUyGXFxcXTp0oURI0bQq1cvdu7cyXvvvUdERATvvvsuX375Jfn5+dJr2NraUl1dLc3y/1olXKH1aLpd2vbt2/n000/JyspCT0+P48eP07FjR+BB8BMZGYmTkxNRUVF4enoCUFdXx/r169m7dy/V1dXIZDL8/PykLEah7TAxMeGVV15hzZo1FBYWSo8fP36cwsJC5s6di6enJ1evXuXo0aPY2dlRVVWFjo4OgwYNori4mN27d0s/Z2BgQHx8vFQXQntvaC3BMogZ5mdOZWUlEyZMYNSoUZiZmaGnp4e+vj6HDx+msrISW1tbRo8ezZ49e7CxsSEgIAAzMzOMjIyk9c62trasXbsWc3NzoqKiWtUIkCAI7Yd2thj+uVbu4fuNNujNyMigsLBQWneqo6ODj48Pu3fv5tSpU0RHR0trKZtqz/cu7bnduHGDK1euEBQUhIWFBQARERGsX7+epUuXUlxczIQJExg2bBg2NjZ89dVXHDhwAA8PD+zt7ZHJZNy5c4fMzEw6d+6Mvb19S57Wf0zbZjQaDbdv32b+/PmcPXsWf39/pkyZQvfu3dmzZw8ymYyoqChpb+SAgACmTJmCq6srPj4+uLq60qdPH3r06EHfvn2xtLRstobzWWBqasratWtpbGyUAiRdXV2cnJzYvHkzwcHBLF26lPDwcJKSknB3d2fVqlVs2bIFc3NzLCwsMDQ05MqVK+Tn59OzZ89n5tq1dTKZjNLSUo4dO8aUKVNITk4mOTmZ6upq0tPT6datGw4ODsCDyZiGhgamTZvGqVOnuHjxIn/729+4ePEif/7zn6UdDJq+ttD6aTQa6Z7n6+vLjh07uHTpEklJSQBcvHiR/Px8OnTowNixY9mwYQMjRozgq6++QqVSAeDm5kZmZiaHDh0iICBAylRt2gZaY8aOCJifIWq1GkNDQ9LT01Gr1YSGhgIPGm9OTg6bN29mwYIFJCUl8d5779GxY0eWL1/OihUrcHR0xMnJSeq8lpaWcvnyZZKSksSNThCEJ6ppJefq6mq+//57bt68ia+vb7P7jfZ52grYc+bMabaMBKBDhw7Y2to+85X9p0yZQkpKCra2tlIFbFtbW9LT0xk9ejTJycnY2Njg4eFBQkICeXl5TJ06lVOnTtG7d29MTExYt24doaGhba6olTZDSluz48MPP+TkyZMMHToUHx8fzMzMsLOzY9GiRfj7++Pp6YlMJkOlUlFQUMDixYvp06ePlLKvnalvWln7WVFfX09eXh65ubl06dJFqhJuaWlJSUkJBw4cYOjQoSgUChobG3FwcCA2NpaSkhIWL15MQkICtra2VFZWEh4eLgVYQuuXk5NDjx49uHnzJv3792fw4MEYGhpiZmZGbm4umZmZ9O3bF3hQMyIqKgonJycaGhq4ffs2Xbt2ZfLkyTg5ObXwmQj/Le33ck1NDfv378fFxYXFixcTHx+PnZ0dJ06cYPv27axatYqYmBgmTZpEXFwccrmcPXv2cOXKFTw9PXF1dcXT05P4+PiWPqV/mwiYnyEymYzq6moOHjyIgYEB4eHhUgCsUqk4ffo0bm5ufPfddzg5OeHp6Un37t25d+8eCxcuZOfOnSQkJGBkZMSBAwewt7cnMjJSzDALgvBEae8nS5Ys4fXXX6empgZvb29UKlWzvVu1z1Mqlfj5+VFeXs7ixYuxsrKSin3Z2toSEBDw9E+ilVCr1ZiZmXHgwAGuXLkizQTo6Ojg5ubG8ePHyc7OxsfHB5VKhUajwdLSkm7duqGrq4uHh4e0/Y+bmxuurq4YGxu35Cn9x1JTU1myZAk+Pj44ODhgbGzMsWPH6NSpk1Sky8vLi4MHD5KVlUVkZKQUHMfExHDv3j169eol1lnyYFlDTk4Op0+fxtnZWRqcUigUuLi48OOPP3L//n1iY2OlAQUzMzMSEhIYMmQIjo6OANjb2zcb2BJaj8ftR669L2RnZ3Pw4EFeeOEFOnToADwYLKmrq2Pfvn2YmJjg6+sr7a7i6+tLfHw8PXv2lNalNk3tFtqejRs3MmzYMO7fv09BQQG5ubnk5eUxYMAAXF1d2bx5M0FBQXzxxRfS0o2ioiIWL15MTU0NUVFRWFtbS+2nrRAB8zNEo9GgVCo5ffo0+/fvZ8CAAVLpfpVKRWlpKWfPnkVPTw9fX18aGxsxNDQkNjYWR0dHbt68SUxMDKamplhZWWFlZYWDg4O48QmC8MTt2rWLuXPnMn78eN5//318fX2bBWoPD9TJ5XLi4+MpLS1l27ZtZGdnEx0d3aYLVD0J2lTkuro6NmzYQOfOnbG2tpZmmb28vPjhhx+wsLAgICAApVIpFYUMDQ1ttleuvb19qw+WH9fZ37t3L6mpqVhYWBAYGEhwcDA7duygsrKSwMBATExMAAgMDGT27NlYWVkREBCAXC5HV1dX2mbrWaf9zDk7O7NmzRrgwTXTFvXSLvOaM2cOvXv3xtzcvNnPKpVK6f0R2wb9565fvy61yScVdD78Otoif/BgO7Vr165RX18v1XgICwtj8eLFdOjQgZCQEOm5lpaWXL9+nQMHDtCnTx+pkjz8MxND+7tEn7FteNxkWHl5ORMmTCApKWjqgIwAACAASURBVIm//e1vxMbG4unpyYoVK7CxsSEsLAxdXV0OHDjAzp07qa2tJSMjg08//RQTExPeeuutNltcUwTM7dTjOg3woPPk4eHBtGnT8PX1xcvLS+ocOTg4cPr0aS5cuEBsbCxGRkbSKKGHhwfJycmYmpoCSMGyIAjC/0J7/9HSdqp27NhBaWkp77//PlVVVVy/fp2bN29y+/ZtKQW76Ze59udiYmIIDAwkOzsbCwsLbGxs2nUHTVts6tfOUSaToVAouHDhAqdPnyYlJUXaZ9nGxobCwkK2bduGj48Pzs7Oj51daivXUEdHh9LSUu7cuSN9XwUHB3PkyBGuXbuGl5cXKpUKU1NT1qxZg6OjI/7+/sCD77Vz586RnZ1Nr169Hun0t5Vr8N96+LP4MG1qu5GRETU1NezevRsDAwOCgoKk/3d0dGTz5s3cvn2b7t27N/tZaJ1rE9uC69ev89Zbb5GRkYGfn580c/e/tkuZTEZhYSGZmZk4OzujUCi4efMm7777LosWLeLMmTPMmzcPV1dX7O3tsbS0lPbY7tKli7Qu1cTEBLVazZYtW9DT05OW/Gl/R9M/hdavaf2QpoqKivj222/54IMPsLe3x8DAAA8PD0pLS1m9ejWDBg0iJCQEX19fzp49S35+PufPn+eNN97gww8/xMTEpM3eS0XA3M40XfunDXa1tA3UxMSEgoIC1q9fz8svvyytQTI2NqahoYHDhw9TWlpKp06dHvnAtNWGLghC69L0XqVWq6V9GS0tLZHJZBw+fJjc3Fy2bt3K5s2bOXz4MLNmzWL79u14e3vj6urabGBQe1/SaDRYW1sTExMjFaxqr5oWnNLOGP8Sbcd2zZo1mJiYSDPHMpmMwMBA9u/fT1JSkvS8plrzNXx4hqy0tJS+fftSX19PWFgYSqUSHR0djIyMSE1NBaBjx4506NCBjIwMLl++LKWjw4M9p5vutazVmq/B/+rhz+L+/fuRyWTNZoibkslkhIWFcf78eTIyMjAxMcHLywtA2ratf//+T/MU2p2mfS2NRoO5uTk9evTg4sWLLFy4EBMTE3x8fP7jdvm4PtyHH37I9u3beeWVV7hy5QqffPIJKpWKf/zjH4waNYp79+7x448/4ujoiKenJ126dGHWrFmo1WoiIiKkPqStrS2hoaH06dPnyVwE4al7+F7w008/UVFRgUqlQi6Xc+/ePQ4ePIiTkxOBgYFoNBoUCgXGxsYsX74cQ0NDIiIisLe3p3fv3iQlJTFo0CBp6UtbjiFEwNyONG2ImzZtYsCAAcTExDx2JrhDhw6sWLGCkpISQkNDpc6Bs7Mzp06dwsbGhvDw8EdGg9tqQxcEoWVlZ2c/dlZk9+7dDB06lJ9//pmNGzdy7do1unfvjq+vLyYmJpSUlJCYmEhMTAyjR4/mwoULHDt2jP79+/9iFg2035kstVotFbDSdmq+/fZbNmzYQH5+Pj4+PlIHVkt7vV1dXVEqlUyePBkvLy9cXFyQy+UYGBgwYMCAxwbLrc2RI0doaGjA3Nyc+vp6KSVUy8DAgMLCQo4ePYqnp6dUoMzNzY0zZ85w8uRJXF1dcXR0xN3dnenTp2NhYUFYWBhyuVwagGjP6yz37t2LRqORKqZrz3Pt2rWMGDGCo0ePsmLFCjIzM5vtlap9rvbaeHt7c/fuXWbMmEFgYCBWVlbo6upKgXZ7voa/lV/LGDE2NqZr167k5uayfft2FAoF/v7+//Z1fjhY0Q6yyeVyfvrpJ/r164e+vj5yuZz/+7//w8zMjA0bNrBo0SLu3LlDY2Mjfn5+mJmZYW5uzrx58wgPD8fFxQX4Z7X0x/0uoW3Qvmfp6ekMHTqUY8eO8dNPP3H69Gl69eqFUqnk2LFjlJSUEBwcLC1nycvLY+PGjWRkZJCSkoKFhQUymUwaxG0P6fgiYG5HZDIZ2dnZHDp0iFOnTjFw4EC6dOnySIcCHqw18vT0ZO7cuejo6ODv749SqUSpVNKpUyc6d+7cbjucgiA8XbNnz2b27NkMGDAAhUKBTCYjKyuLQ4cOkZeXR+/evRk7dix2dnbMmjULZ2dnKa2rV69e+Pn54ebmhrm5OVu2bCE+Pp6wsLCWPq2nZvPmzRw6dIiwsDBp9B8gPz+fjz76iCtXruDm5sa8efOoqqrCx8cHY2PjRzopSqWS8PBw7t27x/79+8nMzGyWNvtLS3lai+PHj/PBBx+go6NDTEwMcrmcY8eOsXTpUm7duiUVhQsPD2f58uXU1tYSGBgobRPl4ODAnDlzMDQ0JCwsDCcnJ3R1denRowfW1tZA+08fzcjIYMGCBcTFxTULhI8cOSKlWr799tvExMSwdu1abty4ga+vL+bm5lIQpL02pqamdOrUicrKSnbs2EF2dnazqrft9Rr+VrQDYdr74549e6ipqcHKygqlUillDfr4+HDnzh1mz57NoEGDpDXk/4pMJuPMmTNMnDiRrl27SgNr5eXlHDt2TKpZ4+7uTn19vTTz/NFHHxEfH8/333+Pu7s7gYGBBAQEsHTpUgwMDOjUqZOYXGknLly4wM6dO8nLyyM5OZnPPvuM0NBQZs2ahY6ODnFxcajVavbt20dBQQGdOnVCo9GQmppKeHg4rq6uREREPJKd0h7agwiY27CHOzclJSV89913LFy4EHNzc8aOHftIWpmWRqPB3d0dIyMj0tLSSE9PJzIyEiMjI/T09P6tdXGCIAi/5O7duzQ0NKCrq0twcDCvvPJKs8G7NWvW8M0331BdXc2bb76JtbU1AQEB3LlzhzVr1pCSkoKRkRHV1dVkZ2dz8eJFPv74Y4qLixk1apQU4DwLVq5cSXR0NM7OzjQ2NlJVVcV7770nZQNNmjSJHj16YG9vz7Zt21AqlYSGhj5y/9YG0B07diQwMJC1a9diYGCAjY0N+vr6rTpYhgcBb05ODhcvXsTe3p60tDQ++eQTTExMWLJkCSUlJXh4eEj7em7cuBFXV1epGuuZM2fIyMigoKAAR0dHvLy8CA8PbxY4tnfOzs70798fKyurZmuWN2zYwM2bN/nTn/6Eubk5zs7OWFtbs2/fPvT19QkJCfnF9hQVFUVsbCxqtRpXV1dRIO2/JJPJKCoq4p133mHhwoXcuHGDZcuWceXKFZKTk6Vg2sjICFdXV/bv309WVhZJSUmPnWV+3Czvxx9/zJ49e6irq8PCwgJra2tMTEzYsmUL+vr6xMTEYGhoyI4dOzhw4ADTpk2jY8eOWFhYsGrVKiorK3FwcMDR0ZH+/fuTmJjY6u8bwuM9boD04MGDjBs3jpKSEsaMGYNKpZIyCBYtWkSXLl2kz/oPP/zApk2bWL58OUePHuX1119n6NChv7iUo60TAXMbpG3kDzd0Q0NDNBoNp06dws/Pj549e/5qqo5MJiMgIIDY2Fh27txJUVER9vb2mJubt/nUCUEQWk51dTU9e/akvLycuLg4KVCeN28eRkZGqFQqnJ2dOXfuHAqFgqFDh0qdu+DgYBYtWoRGoyEmJoZLly4xb948fvzxR6Kiopg9e7ZU9Ku9096/u3btirOzM5WVlejo6GBgYMDevXvZvHkzERERJCQkAODv78+xY8fIzs7Gw8MDW1vbZp3mpunqVlZW9OzZE1dX11bfwWlsbJQq7drb27N7926qqqooKCjgiy++YOTIkbi4uHDgwAGKi4vp0qULwcHB7Nu3jwsXLkh7Kq9evZqkpCQ6d+7cbJ3ls5Q+qj3XtWvXsnnzZsLDw9HV1WXjxo3U1tYyaNAgqY/RoUMHdu7cyb179x4blDVtV8bGxnTo0EEEy/+DpUuXMnLkSCIjI5k4cSJ9+/bFysqKRYsWERERIQ2Y6ejoYGpqiqmpKfPmzZNSYB9ux03/rp2d1tPTIyMjAysrKy5duoS9vb00EJWWlsbIkSNpaGhgzpw5GBgY0L9/f/T09Ni2bRt37tzh1q1bJCQk4OTkhJ6eXrtfvtAe/VIMAeDt7c3Fixe5d+8egwcPRl9fH4CoqChWrFhBcXExcXFxhIeHk5CQgIuLCx4eHkyfPl1a/tle76ciYG4DHm582kb+008/MXnyZE6fPk1hYSGBgYE4ODhw8+ZN0tPTee6556TKhQ833qbFJExMTOjWrRshISHY2dm1y4YuCMJv79q1a1hYWKBUKpHL5Sxbtozu3btLW5KMHTuW7OxsUlJSMDU1RU9Pj0WLFhEXF4e9vT1qtRpDQ0MMDQ2ZO3cuXbt2xd/fH3t7e0aNGkVycvIz1UFreo4ZGRl88sknGBgY4OPjQ1hYGDt37sTa2prIyEipY2NlZcWuXbtobGwkKioKuVz+ix1pAwMD6edaI7VajUajkdYWV1VV4eDgQFFREWvXrkVXV5fRo0cDDzp6169f5/jx49jY2ODm5oanpyeZmZmsX7+epUuXoq+vz3vvvUdwcDDAIynG7UlNTQ0KhYK6urpm56j9c+PGjRw7dgyVSiXNwM+fP59+/fphYWEhrW+9cuUKR48eZejQob86+C785w4fPszly5dxd3enuLiYlStXUlNTw+LFizEzM8PIyIjS0lK2bt2Km5sbUVFRzYIcpVLJpUuXqK6uJioq6pH34fjx44wfP54uXbpIyxLKy8u5e/cuwcHB1NfXs2jRIgYNGoSxsTGrVq0iJCQEV1dXbt68yaZNm6ioqCAzM5OlS5cybtw4/vjHP+Lp6Qm0/+UL7ZW2DW3atIkJEyZw6tQpysrK8PX1RUdHB0dHRxYsWIC/vz8dOnSQAmwnJydmzJhBQEAA7u7uqFQq/Pz8CA8PRyaTNStC2R6JgLmVUqvVlJSU8PHHH+Pi4iLNqGg0GsrLy3n//fdZv349Xbp0QSaT8f3331NTUyOlzpw8eZL8/Hy6du36L7cbAdDT08PAwKDdNnRBEH5bmzZt4ve//z29e/fGwsKCoKAgNm/eLBXxksvl+Pn5MXPmTPz8/PDw8MDBwYHz58+zZ88eXnrpJen+4+/vz/Lly9FoNMTFxWFnZ4ehoWG7KBzyr2iDxMcNcm7btk1al2tjY0NFRQXp6el4eXnh5uYGPEhbzs3NJT09HTs7Ozw8PNrs9dK+12fPnuXzzz8nPz+f6OhofH19OXLkCLW1tXTv3l0KBiwtLTl+/Dh5eXl069YNBwcHunbtSmRkJP379+f1119vVhCtrV6XX6LtH/zlL3/h2rVrREZGIpfLkclk3Lp1i6qqKoyMjIAHn7G0tDSKioqktYfHjh1j7969Uq2BxsZGFi9eLKVcC09OTU0N8+bNw8zMjICAAIyMjDA0NOT06dPcuXOHjh07cubMGb755hvu3LnDa6+9hqurK/DP7b+MjY05cuQIenp6xMbGPnLfKC4uZtq0aVy/fh0nJydsbGyora1l7969REVFMXjwYFatWkVubi4WFhbcuXOHiooKaQaxtLSUCxcucOnSJT788ENiY2PR09N7ZgYs27Km71FtbW2z5VCFhYW8++67bN26lZiYGGpra5k+fbpUCNHJyYn8/Hw2bdokLY0C8PDwYMuWLVRVVT12b/r2npovAuZWSJt6ZmRkxHvvvYdSqZRK92s7D6mpqcydO1dKL1u1ahUFBQXExcXh5eXF3bt32bVrF0FBQdjZ2bX6Yi6CILRtLi4u7Nixg+LiYhISElAoFDg6OjJlyhQiIyNxdnbGycmJ7Oxsdu7cSZ8+fTAxMcHR0ZH58+fj7OyMj48PjY2NKBQKUlJSSE5ObvY72nsnrbGxUQpwCgoKqKyslKrWmpqaUldXR3p6OoaGhgQFBREVFcW6desoKysjKCgIY2Nj4ME61dzcXHr16iVVQm6tfm3tZWNjI5MmTeLLL78kJiYGf39/rK2tpfTTzMxM9PX1pRljlUrF/fv3ycjIoLGxkZCQEHR1dXF0dJTS06H9tiOZTIaenh7Lly+nsLBQ2q/3//7v/5g+fTpbtmyhsLAQFxcX7O3tqa2tZf/+/RgYGBAZGYmvry+LFi1i586dXLt2jUWLFpGdnc0f//jHx+62Ifz3FAoFCxcuxN7envDwcOBB+y0uLmbbtm2kp6czb948/P39CQoK4vjx41LhP20AI5fLuXbtGnv27GHQoEGPtGs7OzvCw8PZsmULJ06coGfPntja2pKamkpubi7PPfccYWFh/PzzzxQUFJCbm4ulpSUhISEYGBgQFxdHYmIiw4cPx93dvV1nZLQXTd+j8vJyZs+ezcmTJ1GpVFhaWlJXV8fu3bspLy9n2rRpJCYm0q1bNzZs2EB2draU0RUcHMz8+fOlbaK09+m+ffuSkpLyTC69EAFzK7NixQpSU1Px8/NDX18fJycnZs6cKRV8AZg1axampqYMHDiQadOm8fbbb+Pj48P48ePx8vKSOlfnz58nLS2NgQMHimBZEITflFKpxNramqlTp9KpUyccHBxwc3MjKyuLtLQ0nnvuOXR1dQkMDGTWrFmYmZkRGhqKtbU1ZWVlTJs2jddeew2lUgkgBX/P0myGjo4O5eXlfPzxx0ybNo1du3Zx6NAhvLy8UKlU+Pr6sn//fq5fv06HDh1QqVSYmJiwfv167Ozs8Pf3B8Dc3JzExMRWHSw/3PnesWMHt27dkgrMyGQyrl69yqJFi/j0008ZMWIE3t7eUkVgPz8/fv75Z65evSpdC3gQJJw6dUqq5ttUe+/sa2ceHR0d2blzJw0NDRw8eJCysjI+/fRTDAwMOHz4MNu2beN3v/sdAQEBHDx4kJycHLy9vfH39ycqKgqNRkNOTg4uLi7MmTNHek+EJ0N7TyssLCQrK0taaqKrq4uBgQFnz57l4sWLbNiwgUGDBpGYmIidnR1bt25l48aNVFdXExgYiI6ODrW1teTm5hIXF/fY5RVOTk5YWlqSlpbGgQMH6N27N3Z2dixYsIB+/frh4uKClZUV169fZ+/eveTn5ze7D2szMrQVvIXWTXt/mzVrFu+88460naCNjQ12dnbI5XJqa2uJi4vDwcGB+fPnM2bMGPz8/Lh06RINDQ2EhoaiUqmQyWRMnjyZlJQUaVmVUql8ppZFNSUC5lbmyJEjzJgxg65du+Lg4ICvry+7d+/m0qVLREdHY2xszLVr19ixYwfLly8nLy+PTz75hPfeew+VSsX169dpbGzE2dmZsrIy1Go10dHRUiMXBEF4EtRqNdB8ts7Ly4tjx46RkZHBc889h1KpxM/Pj+nTp2NnZ0dgYCBmZmbU1taybNkykpKSMDc3x8nJibKyMjp16vTIFint+b71cKfjzJkzvPnmmxgYGDBu3DhSUlLYtWsX58+fx9vbG2tra/T19dmzZ49U7drb25utW7eSk5NDXFycNPv0uNdvTbTHtWvXLt58803OnDlDbm4uCQkJ0u4O69at4+DBg/z1r3+Vfk6tVksz8WZmZuzZsweNRkPHjh2BBwMtCQkJ0qzds0KtVkuzPg4ODly6dInjx49z9epVxo4dS1hYGNHR0Xh7e/PDDz+gp6dHWFgYxsbG7Nu3D7VaTceOHbG3tycmJoZevXrRvXt3KTVbBEtPjrbtnzp1itzcXMLDwzEzMwPA1taW8vJyLl26RGhoqDRY4eDgQK9evbh79y5TpkzB0dERf39/8vPzpaUrD9MOSmkHj2bPnk1BQQF2dnZoNBoaGhrw9fXFwcGBwMBAUlNTsbGxoW/fvlJG48PHLLRuarWaqVOnsm/fPj7//HPGjBlDWFgYjo6O0nPs7OywsLBg9uzZbN++nT/96U+MHTsWgB9//FFaoxwUFERxcTFdu3aVlr08y+vWRcDcyoSHh7Njxw6uX78ubfPk5+fHpEmT8PX1xdfXl5KSEjIyMvD09GTFihVSAYaKigomT56Mjo4O3t7e+Pn50atXr0dufIIgCP+LpqnD2uJCWt7e3kybNg03Nzd8fX2xtLTk3r17rFy5kpSUFIyNjQkODmbx4sXcuHGDpKQkLC0tSUpK+rf3E23r1Gp1swBHW8E2Pz8fhULBX/7yF+zs7Lh+/TrLli3j9u3bKJVKoqOj8fT05OzZs1y6dAkrKytcXV2JjIwkMTERe3v7Zr+nNd/36+rq+Oabb1i0aBHDhw/n3XffJTY2VpopBjhx4gRFRUV07dpV2ldaW91Vu4VRVlYW6enpUiohPJgVa6+VWn+Jdp3y2bNncXZ2xt3dne3bt5Ofn88777wjzT6qVCoqKipYs2YNI0aMwM3NjTNnznDgwIFm11BbLA7a/9rEp03bNk1MTJg6dSqRkZF4eHgAD95HU1NTrly5ws8//0z//v2BB9kDBgYGODs7k5SURI8ePWhoaMDV1ZXg4GDpHvy4qvjwIEhycnIiLS2NkydPolAocHBwICAggIaGBoyMjOjXrx/Dhg2Tql8Lbc/du3eZNm0aL774Iv369UMulz9S+FEmk3H79m1mzJhBYmIigwcPBh7UITlx4gS3b98mNjYWMzMzevToIQXLzzoRMLcC2oas7TR5enoyadIkgoODcXNzw8HBgStXrrBr1y569uyJh4cHV69e5fTp03h5eQFw+/Ztxo8fz7Vr10hJScHBweGxN1BBEIR/1y/dO7SpgF999RVLly4lPT0dW1tbLC0tsbe359atW6xfv56UlBQMDQ0JDg5m+fLl1NXV0alTJ3R1dbGzs8PMzIyQkJB/+fvaE22NCh0dHXJycpg8eTL19fW4uLhgbm5OaGgotbW1fPjhhyxcuJCRI0eiUCg4efIkXl5eODg4YGNjw+rVq7G2tiY8PBwLCwvMzMxa9Yzyw7Kzs1m+fDlffvklffr0wdjYGFNT02bPKSgoYP/+/Tg6OuLr6yudW2FhIUePHsXDwwM3Nzfc3d2lrbW02sp1+G89/F7fu3ePzz77jPT0dLp37y6tUb5w4QKmpqYEBQUBDz67d+7c4ezZs8TExGBpaYm7uztubm6PvYbt/Tr+lpruc92UNqVVpVJx+vRpjhw5QlJSkpT+bGFhQX19PYcOHUKj0RAcHIxGo6GxsZFOnTphbGxMx44dUSqV1NTUsH79egoLC7G3t29W1K7p7wOkJQ2nTp0iIyOD+vp6KagCpMmVtnQfEZrLzs5m3bp1vPDCC7i5udHQ0MD27dvJyMhg69at3Lp1i4CAAJRKJRMmTCA2NhZ/f38KCws5efIkw4YNw9nZuVmRv2fhe/nfIQLmFtTY2NjsC0l703JycuLMmTMcPHiQuLg4zMzMiI6OZurUqRgaGhIfH4+npydFRUVMnz6dQ4cOsWzZMlxcXJgxY4Y0UqklGrogCP+Jh4sjPdzxu3jxIr///e+pq6sjISGBrKws1q9fT0NDA5GRkYSGhrJgwQLkcjnR0dHo6+tjZGTE5MmT6d69OzY2Nnh7ezcLlpv+vvZMW8xqwoQJfPbZZ7i4uODr64u1tTXm5ubo6+uzYMECcnNzmTRpEomJicjlctauXQtATEwMjo6OBAUF0bdv32bbeLSl63f06FE2bdrEu+++K2UWXLhwgcLCQi5fvoyzszO+vr6kp6eTmZmJlZUVDg4OVFVVMXv2bNLS0oiPj8fR0RFvb+8WPpvfzi91Vh9+TF9fn5qaGs6fP091dTURERF4enpy5MgR8vLyCA4Olvbb3rdvH1evXuW1115DLpdjYWHRrq/h06bt22n7dFlZWZibmz9SKEkmk+Hj48PkyZOxsrLCx8dHWjusUqm4cOECu3btYuDAgSgUChQKBbq6uqxcuZJOnTqRnZ3N8OHDuXDhAtu2bePEiROEhYVJ6d1NaduRh4cH7u7ubNiwgYiICLp37y4dS1u8jzyLsrKyMDY2RqlUPnJ/sLOzY/PmzaSlpbF//34mTZpETk4O2dnZnD9/nk2bNqGjo9OsOva+ffuYPXs2Pj4+vPHGG8/k9/K/QwTMLUQ7yyCTyTh69CiHDh3i7t27WFpaoqurS2RkJJMmTcLe3h4fHx9MTU1pbGxk2bJldO7cmQ4dOtCjRw+SkpIIDQ1l2LBhDB06FH19fTE6KAjCf02b9iqTycjIyGD+/PkcP34cuVyOvb09Ojo6bNq0iWvXrrF06VKioqIYOHAgubm5HDt2DAcHB7y9vVEoFCxYsIAePXpgaWlJQEAAxcXFdO7c+ZGZxPbscffjo0ePsmTJEiZNmsTIkSNxd3eXgsby8nLmz5+Pq6srL774IgAbNmygsLCQ8vJyKW1WmzrbVu/3+fn5ZGVlkZubS2VlJd9++y1paWksX76cVatWcfjwYYKDg+nevTvnz59n6tSpHDlyhFmzZnH//n3Gjx/f7otRNS20VFVVJQVTAFevXmXy5MlSmi4gpaifOXMGf39/HBwcUCqVbNy4kZ9//pmGhgbOnj3L3Llz6d27N506dRKzyE+QRqORllrIZDLS09N544032LFjBz179mwWyMpkMjQaDVZWVtKAmLGxMQEBAQAYGhpib2/Piy++iLW1NfX19SgUCiIiIli5ciWVlZWUlJTQp08fxo0bR5cuXViwYAGNjY0EBASgr6//i+nZTk5ODBkyhOeff168/23M2rVrGTNmDDExMTg7Ozd777QD27GxsRgZGVFUVMTAgQPp06cPI0aMYOTIkdy/f58NGzYwcuRIoqOjCQoKQqVS8fbbb/O73/1OtIVfIQLmFqC9iZWWlvLOO+/w/fff09DQwMKFCykoKMDJyQk3NzcqKipYvXo1cXFx2NjYEB0dzcqVKykuLiYyMhIDAwMsLCxwdnbG0tJSKsIj1hsJgvDfkslkFBcX8+6770oBcVVVFV5eXlKxmGXLlmFhYUFKSooUsNnY2HDo0CFqamro3LkzwcHBbNiwgTNnzkjbUHTv3v2ZCZabdp4BysrKpKB40qRJ6Ojo8Oabb0pBkfY66unpSXtd6unpcfb/2TvzgJyy/oF/etqUNq2UotJCkVQqWlBiyD7WsY0xCMMM4/UyY6zDDLKMQbaxzRjrRIUokp2ikUKWkZpQpFV73d8ffs99i8xqq+7nH3qee89z9FWG7QAAIABJREFU71m+53zP+S5Xr3Lq1CkmTpzIsGHDXjgJrGkLHPn8Z2hoiJKSEidOnODIkSO0bNkSb29vhgwZQv/+/Vm/fj2Kior4+/vj5+eHk5MTzZo1o0uXLsyYMQN9ff1aayooD7KloKBAamoqCxcuJCwsjBs3bqCpqYmRkRHJycls2bIFmUyGs7OzePqooqJCbGwsjx49wsvLCysrK5KTkzl//jwAjx8/ZsCAAYwePbqKdYLEv0dBQQGZTEZKSgoTJ05k27ZttGnThpSUFIYMGSJG/n/+HhcXF1JTUzlx4gS5ubk4OjoCz5TmTZs2YW1tjba2NmlpadSvXx8bGxsWLlzIw4cPmTBhArq6uhgYGKCgoEBYWBjm5uYvzb8u/0xNTa3awI0S7zZ2dnaEhYWJ/aSyf7F87S/PQNGtWzfs7e0xMTERYxjIMwt07doVTU1NmjRpQps2bTA0NJT6w5+g9OeXSLwK9u3bh4WFBY6OjuIkv2PHDlHAGRkZcfbsWT799FOKi4sJDAzkv//9L/v372fv3r2YmJigo6PDpEmT+PLLL5k4ceILKUMkRVlCQuLfkpmZyezZs1FSUmLXrl1VTvHksktLS4u4uDjxM4AWLVqgo6PDgwcPUFBQQElJiZkzZ3L9+vUqQcFqo5JT3TvJTTJTU1NZunQpysrKjB8/HgsLCzEQS2WzzcryOyAggDVr1rBw4UIqKiqYMmUKHTt2fOlvvYvIA8M9j/xkTUNDg6FDh9K9e3e0tbVfmL/Mzc1JTk4W37eyTx3UnHr4J8iD9MydO5c9e/bQtWtXmjRpQlhYGLGxsWzcuBFHR0e6d+/O6dOniY6OpkuXLgB4enpy6NAhTp06xfnz53Fzc6N79+4kJSUxatQonJ2dxd+pzXX4png+3VJwcDAzZsygR48efP/990RHR5Oeno6RkdEL18r9hWUyGZMnT+bcuXMsWLCArKwshg0bJm50FBYWoqyszO7duwkKCsLDwwMfHx+uXr0qKjkAH330EYcOHeLw4cM0b94cExOTP2xjac1YcygvLweeyYapU6cyffp0vLy8xJRkf4WUlBQSEhLo16/fCwEiQeoPf4ZUO2+A+Ph4goKCCA0NFYVjXl4ep0+fplevXhgZGbFnzx7+85//YGNjwwcffCAuQv/73/8SHBzMlStXEASB3r17c/r06VpviiYhIfFmkcucU6dOcfnyZT7++OMqcqaiokK8ZvDgwdy9e5fQ0NAqSlF5ebl4DYC3tzfjxo2r8ju1aYH+vK/385/v3LmTnj17oqysTNeuXREEAUEQsLa25tGjR1y6dAn432IoPT2d8+fP4+zszLp16wgKCuLkyZP07t1bLPddrz/5O8r7RWxsLDk5OeJ3ULW+tLS0Xlio3bx5E4BOnTq9tG7f9Xr4N9y/f5+ePXty+PBh9u3bR2BgINOnT2fQoEFcvXqV6OhoAIYMGYKqqirHjh0jPT1dvN/Kyor09HS2bNlCSUkJrq6u/PDDD6KyXBfq8E0h77u3b99GEAQsLCw4cOAAS5YsQVtbm6SkJFE5qU4hkX+mrq6Or68v3333Hba2thQWFtK2bVuaNGnCnj17iI2NJSQkREwfNX36dDIyMjh9+jSlpaVieePHj+fy5ctERUXVCHkh8deQb7KGhobSoEEDtLS0OHz4MPfu3fvD+y5cuMDFixdZvHgxffv2pVGjRowcOfLNPHQtQ1KYXyGVd/oq06pVK7p168b169c5fvw48MwfKS0tjZycHN5//31Wr17N5MmT2bp1K23atCErKwuA3r17o6WlxenTpykrK0MQBHR0dF76WxISEhL/hMp5cVu2bPlCHlt5ZGcAW1tbhg4dyty5c9m+fTt3795l//79JCYm0rNnzzf+7G8LeZ2FhoYyb948MjIyxM+zs7PZv38/U6dOZenSpfj6+mJpaSnmT9bV1eX7778XTbeLi4tZv349x44dIzc3F2VlZWxtbYFnvmmVf+9dRu4TGRsbS58+fZg7dy5Hjx4Vv3seuWL95MkTHjx4wP79+/nkk09o2LCheKr+fPm1HWNjY0xMTHB0dKyiZJWWllJWVoahoSEABgYG+Pv7c/v2bcLCwoBnfeXBgwe0b9+ehg0bkpubS0VFBfXr15dMLl8BlTcE5WzZsoUPP/yQCxcu4ODggI2NjbgJFh8fL6b+lH/2MsrLy2nVqhXdu3enSZMmpKSkoKysLEbHlwd0lUfVHzRoEJs2baqiNPn4+GBiYsLt27cpLi5+Va8t8YaRzwtyMjIyGDJkCIsXL+b06dOoqakRGRnJyZMnKSkpeWk5mzZtYvny5SQmJrJ27VoCAwPR0tKqth9L/DGSSfYrRCaTkZWVRUZGBjY2NsCzyUtJSYnBgwdz7do1wsPDadOmDUZGRrRs2ZK5c+cyduxYPvroI9G3786dO4SHh+Pj44OtrS0HDhyQzK8lJCReK4IgUFpaioKCAjk5OS+YD0JV88Pp06dTVFTE9u3b+emnn8jKyuLTTz/Fx8fnbTz+WyM6Oppp06bRuHFjrKysxJyWV65cITk5mVatWonXyuuvdevWDBo0iMDAQDw9PbGzs+PWrVvUr1+fRYsWveDnXdmkvSawe/duli1bxvvvv8+QIUNERaG6PgVw/vx59u/fT0pKCvfv32f06NEMHTr0TT/2O4HclP3DDz/km2++ITo6Gm1tbebPn09kZCRGRkakp6eTn5+PhoYGAwYM4NatW2zcuJGEhATu379P/fr1mTVrFubm5lXKltYN/57qNht69+7Nzp07OXbsGNbW1ujq6iKTySgoKCAjI0NM6aWoqEhKSgqlpaWiEg3/GxeKiopkZWWRmJiIlZUVZmZmLFu2jJCQEL7//nt++eUX+vbtKyo7s2bNwtXVleDgYCZMmCD6s65duxZNTc03UBsSrwO5ZYCCggJ5eXloamoSGxtLeno6O3bswNTUlPHjxzNz5ky2bduGi4sLzZs3r1KGXI7MmzePnJwcUSeRrEv+OZL0fMXMnj2bgIAA8vPzqaioEBc6DRs2pLS0lMuXL4u77QMHDkRRURFzc3PxuvT0dFavXk1iYqIo8OSRFaVTZQkJideFgoICKioqNGjQgMzMTNFcuPJOtHzBffjwYZ4+fcoXX3zBrl27mD9/PmfPnhWVxdq4e/28/JX/bWFhgZKSEtra2ly8eJEbN24Az06BCgsL0dXVBf6XGUHOe++9x7Zt2/jkk0+wtrZm0qRJhIaGYm9v/4be6N/zshOz48ePM3jwYD7//HOMjY3FKM4vU9isra3x9fVlyJAhHDt2TFSW6+KcJz9xd3V1pU2bNuzatYvOnTujqqrK2rVrmT59OkFBQQQEBHDw4EFkMhmfffYZU6dOBcDJyYn169eLynJtHItvmuf7+fr164mIiBD7p46ODh988AFRUVHExMQAiIFdi4qKaNasGbm5uUyfPp2ePXty9+7dKuXJx8WKFSvo2rUrgYGB9OnTh2+//Zb8/Hw8PDxwdHTkxx9/pKCgABUVFUpKSlBUVGTKlCls2rSJxMREsTx5cLG6OH5qA3Jl9rvvvsPd3Z3c3FzOnTuHgYEBpqamlJWVoaKiwrx588jMzCQkJIT8/Hzgf+NdLkcaNmwoKsvywJKSsvzPkKJkv2JsbGzYunUrurq6Yi6zAwcOEBAQQEVFBTo6Oty/f5+WLVvSunVrioqK2LRpE6GhoSQmJjJ//nx0dHSYO3cuxsbGwP8Gj9TJJSQk/g1/5NMmn0wbNmzIxo0bUVNTw8nJCVVV1SrXRUREsH37dnx9fdHQ0EBNTQ0TExNkMlmV6L61Dfk7paamoq2tLeZTBnjw4AEGBgbcu3ePsrIyXF1dady4MTt37gSe5U6W319SUsKePXswMjKiYcOG2Nvb065dO9H8Wl6H7zLP+yk/efIEFRUVZDIZT58+JTw8nMePH6Orq8uuXbsICQkhMDCQx48fY29vj4qKSpXy1NTUsLCwwNraukrE8NrYj/4K8ve3sLAgKioKa2trFi1ahI2NDVZWVnh6elJUVMSyZcuIi4vDxcUFd3d3OnXqhJeXl1SHrwi5vJSPx/DwcIyNjQkKCiI6OpoOHTqIBxsODg4cOnSI+/fvY29vj7a2NvHx8URHR6OqqsqECRNQUVEhKCioSuA1Odu2bSM0NJTFixczZcoUXFxc+OKLL5DJZHTs2BFFRUUuXLhAdnY2rq6uKCoqcvfuXby8vMjOzqZLly5iFH5pzVizOX/+PPHx8SQkJDB27Fjs7Oz4/fffiYyMZOTIkSgpKVFcXIyamhqZmZkcOXIEOzs7zMzM/rDNpf7w75AU5ldMgwYNyMvLY/fu3ZibmzN9+nSOHDnCRx99xDfffIO+vj4nTpyguLgYV1dXcaFkaGiIIAiMHj2aTz75BE1NzRqbX1NCQuLdozpF7Pk8nRUVFRgaGpKenk54eDgFBQW0b9+eoqIiSkpKiI6OZuPGjXh6etKuXbsXoiC/64revyEuLo6BAwcSERGBoaEhlpaWyGQyiouL2bFjB8OHD6ewsJBLly5hampKkyZNkMlkrFixgiZNmqCmpoaioiLbtm3j8OHD2Nvbi5ui8OLi/F0iLy8PVVXVKkqYgoICiYmJTJkyhaNHj7J//3709PSwsrJCS0uLyMhIQkNDMTY2RltbGzc3N1asWEH79u1p3LjxH/5ebZ735L6Jf9TO8kji2tra5ObmcvXqVTQ1NcVNFU1NTVxcXDAwMODOnTs4Ojqir6+PkpJSFXNOib9Odest+d85OTksXryYoKAg3N3d6du3LytXrsTIyIgWLVqIctDQ0JBVq1ZhYWGBnZ0dN2/eZPfu3SQnJzN//nymTp2Kjo5OlbLLy8vFzY+hQ4fi4+NDfHw8q1atIi8vj4EDB2JpaYmBgQHZ2dls3ryZkpISpk2bxu+//07Hjh3p2LGjqCxL1Cye38S+ffs2CxYsYM+ePTg5OYkBukpLS7l48SLZ2dm0bdsWmUyGIAiEhIQQFxdHbm4u7dq1k/rBa6RmOUbVEMaNG0dkZCRjx45l2LBhrFmzBj09PQA6dOjA8ePHuXz5spjyoX379rRv375KGc+b70lISEj8Eyr7x+Xn57Nz5040NDTw8/MTzYWfZ/r06ZSXl7N582YiIyMxNzcnNzeXGzdu8Mknn/Dhhx++4bd4+zx48ID8/HzU1NT46quvKCgowMPDAz09PXR0dIiOjmbs2LGMGTOGiIgI7O3tGTlyJLdu3WLJkiXiQqa0tJR58+bh5ORUpfx3VcHZs2cPYWFhBAYGoq+vL8blOHLkCPPnz8fHx4eePXty6NAhli9fTnJyMsOHD8fOzg59fX0KCwvFd9+2bRuPHj16y2/09pAvjhUVFXn69CmXLl2iadOmf5j1YvDgwZw7d46oqCicnJwwNTUVN7969epFnz59qlz/rvajdxW5fJTJZBQWFnLjxg309fUxMjJCRUWFzZs3k5KSQm5uLsHBweJmz9ChQ/nhhx9wc3MTNzJcXV3R1tbm0KFDODg44ODggKenJ/7+/uLJvzxWRElJCRoaGigqKqKurk5WVhYqKiosXLiQvXv30qtXLxYtWoS+vj7l5eVoaGgwaNAglJWVOXv2LB988AFjx44V30OKiF2zkMvR59vM0tISf39/VqxYUSUGgY2NDV27dmXTpk3Y29tjZ2dHbm4uGhoazJkzB5lM9tL5XOLVIJ0w/03+yqmvqqoqurq6REZGMm/ePNF/Sz7JGRgYEBwcTH5+Pt7e3qIQle8qS7vDEhISrwq5LDl9+jTDhw/nwYMHHD16lJiYGMzNzavkY5SfMquqquLh4YGrqyuqqqo0aNAAKysrVq5cSdu2bYG6t0CzsrLi5s2bYgTjtLQ0Dh06hL+/P9nZ2dy/fx9/f38eP37MuXPn0NHRwdraGk9PT7y9vbGxsaFdu3bMnz+/RqUFvHr1KnFxcZSVleHk5CRu5K5Zs4a2bdsyc+ZMjI2NuXr1KuHh4Zibm9OuXTvRj7KoqIiCggK+/vprCgoKGDVqVJ0NSCQfLxs2bCAgIIDLly+zdetWkpOTadasmXj6KL9WPhZVVFQ4fvw4JSUltG3bVjzRlJcnWaP9c+T1tm7dOr788kuSkpI4duwYjo6O6OrqsmHDBo4cOYKtrS19+/YVfYc9PDwICgqivLyc5s2bU79+fU6ePMnDhw+5ceMG7u7u2NnZsXv3bs6dO4evry+ampoEBQXxxRdfEBUVRWZmJo0bN6ZevXpcvHiRoKAg1NXVWblyJf369UNdXZ3U1FSCg4OxtLREV1cXFxcXevXqhYuLCyD5pdZU5HJ0586dHD16lPz8fNTV1dHU1ERLS4vk5GQuXrzIwIEDgWd6RbNmzcjJyWHVqlUcPnyYDRs20KpVK8aPH1+jYl/UVCSF+W8gjzr3PNUtHK2trYmOjubXX3+lS5cuKCoqigPE0NAQTU1NevXqJQb0knxOJCQkXgdXrlxh7969JCYm0q9fP2bPno23tzcHDx4kIyMDBwcH6tev/8JmnZKSEo0bN8bd3Z127drh7OyMsrJyrfZTfhnyRWmDBg04dOgQVlZWDB8+nJ07d/L777+TkJBAQUEB/v7+NGvWjOjoaNLT02nRogW6urro6enRrFkzmjVrBjw7XXjXLYhKS0tRVFTE0tKSmzdvcvnyZfHUODU1lT179jBx4kSuXr3K6NGjuXPnDvPnz2fEiBHiJnB4eDh79+7l66+/prS0lIULF9KkSZO3/WpvjOoi0h4+fJj169ezaNEiRo0aRZs2bdi2bRtZWVnY29tXOxatrKw4ceIE2traODk5vRA1vS6NxVeFvI7v3bvHhAkTuHTpEjNmzGDgwIH4+vqKm1pWVlacPXsWmUxG9+7dUVZWFpXmRo0asXnzZuLi4khJSWHdunVMmzaNSZMmif74HTp0YPny5RgZGREbG8uRI0cYOXIkqqqqHDx4kCtXrtClSxcKCgq4e/cu/fv3x9vbG3gmJ37++WfOnDmDu7u7uKFS+ZBFavuaycWLFxk5ciRXr16loqKCY8eOcf78eXr16kWDBg2oqKjg/PnzlJWV0bp1ayoqKlBXV6dDhw506NABc3NzJk+eTO/evaU+8IaQFOa/gHyxJJPJuHPnDitWrCA+Pp769etjZGT00s5qbm7OypUradGihZg/T16Wra2tmAtN6uwSEhKvgur8lBMSEpgzZw6PHj1i4sSJaGtro6enR0VFBadOnUJNTY2WLVv+ZTn0rit6rwN53cjzm16+fJm2bdsyaNAgEhISOHz4MLdv36Z///4YGBhQUlJCaGhoFdlfmZpQh/LNYUVFRfT19Tl58iQZGRl4e3ujra3Nd999x/79+8WI2F9//TU2NjYUFRURHR2NoaEhysrKZGZmMmDAAD7//HN0dXXrzGloWVkZioqKKCgoUFxcjJKSEhUVFfz4448ATJgwAS0tLSwsLFBVVeXkyZM0aNAAW1vbKvUjH9NeXl507NixxqUYe1eR1/HmzZvJzMwkMDAQR0dHNDU1q8SQ0dPT4/HjxyQlJaGnpyfmUldQUMDa2hodHR0eP35MYmIiY8eOxcfHR3RBkOe/LiwsZMuWLRQUFPDJJ5/w3nvv4eHhgb6+PqdOnSI3N5cPPviA5ORkNmzYwJ07d7h58yZff/01V69eZdq0aWJqquefX6Lmce/ePZYvX06nTp1YuXIl/v7+aGhosGnTJurVq0ebNm3Q0dHh4cOHREZG0qVLF+rXry9uYurr62NrayvKU5D6w5ugVijMiYmJGBoavrby5abSFy5c4OOPP0ZXV5fTp09z8eJF9PT0sLCwqHahamxszG+//caPP/7I+++/T7169V4aVKK287rbSOLfI7XRu8/L2kguf+TRceUTKzxLe5SWlkZaWhpeXl40bNgQgBYtWhAdHc29e/eqTL51RSY9z59Fp64cuTgiIoL79+/j6+tLx44dRX/SytFvu3Xrhpub25t49NeCIAjMnDmTuXPnMmXKFB48eMClS5cwMDDA3NxcPBXZsmWLePIGz1JKHTt2DFtbW5o1a4aTkxNNmzYVy3xXNgtet7yTyWSUlJSwdOlSQkNDadeuHaqqqvz000/o6enRtWtXcZy2aNGC3bt3o6SkhJeXV5VxKK8vebT6urTJ/rrbKDU1lS+//JI+ffrQqVOnKt/JTeLlY/7YsWNkZmbSunVrNDQ0qrSdr68vAwYMoFmzZmJ7VQ6O165dO3bv3s3169f58MMP0dfXB56l/MnIyODkyZN07dqV9957Dw0NDbKysrh37x4eHh6sWrVKdOt7F5HWDS/nZXOKjo4OOTk5DB8+nPLycr755hvWrl2LhYUFhw8fZtCgQRgYGKCoqEhMTAxJSUn4+vpWa+EqWRm8Od6NmetfUlRU9FrLj42N5dNPP2XPnj18+eWXfP/996xZs4amTZuyfv16iouLUVRUrDbf4ZQpU/Dz86vzketedxtJ/HukNnr3eVkbySfSbdu28f777zNx4kTWrVvHkydPABg9ejTl5eWcOXOGp0+fivcMGTKEjIwMgoODgZpx8vm6UFRUpKysjPT0dEpKSoCq+VflZpAmJiZ07dqV69evEx4eDsCYMWPo3bs38EyhKSoqws7OTvy7pnH8+HEOHz5MWVkZS5YsQVVVlW7duqGtrc2BAwcoLS2lb9++mJmZsXTpUnbs2EFiYiLffvstc+bMqZJ7Gao3TX7bvGp593zO28jISDw9Pbl+/To2NjY8ePAAAF9fX44cOSLm05UrXubm5iQlJQF/PA7fpTp83bzuOam4uJj8/HxatWoFPLMKeB555oAePXpw69Ytjh07BiBuEEFVE2mAEydOEBcXVyXA3dSpUxEEgd9++038HXnk87y8PEpLSwEYNmwYCxcuZNWqVYwfP158hncVad3wIvL2ks/LYWFhHDhwgOTkZPGawYMHU1xczJgxY7hx4wY//PADgYGBqKurs2TJEgAcHR3x8PDAwMCAioqKGjmX1Cbq7uqoGuTpHp4nMzOTW7duERsbi4eHBwC2trZ07dqVwsJCNmzYIN7/PCYmJsyePfuFvJMSEhISr4qUlBRGjBjBTz/9RM+ePXF2diY4OJj169eTnZ2NpaUl3bp1Izw8nPj4ePE+Ly8vzMzMKC4upqCg4C2+wduhsszetm0b3t7efPzxx0yePFn0U6yOwYMHY2RkxKlTp7hz5w7AS03j3mUFR55PuTKXL19myZIlzJo1i+bNm+Pu7g48i9Lq5eVFSkoK+/fvR1dXl7Vr16Kjo8PmzZv58ssvuXjxIqtXr2bSpElVFIp3uQ5eBZXTRJWXl1NWVsbu3bsZNGgQW7Zs4eOPP8bKygp4lpPbzMyMmTNnAs8Ur/z8fNLS0l445ZT4d9y8eZMbN24A1a/PsrOz0dXV5fLlywAvmLvLrXZyc3MZMGAA2trahISEkJqaWuU6+SlfbGws3bp1Y86cOcyYMYMhQ4aIZXft2hV7e3t+/vln0tLSqtyfn5//wrMpKyuLMqUub2TWROTt9fvvvzNo0CAWLFjA2rVrGTx4MLGxscCzvhYREUFmZiZff/01Dg4OKCkpoaioSHBwMDExMWhoaDBu3DimTJlS5+KGvItIo/D/KS8vF/2U8/LyquwMdunSBR8fH4qKisQdYHiWQqB9+/YcPHiQlJQUZDJZlROJykg7QxISEq+C6k5B7ty5Q6NGjdi9ezcjR47kww8/REFBgWPHjomnoJ988gllZWVERERUkW8LFizgyy+/RF1d/Y29w7uCgoICSUlJ/Pbbbxw5coSZM2fy/vvvc+3aNWbMmPHCBmrlyMU9e/bkypUrJCQkADVvUVtWViYu9CvPTzY2NnTv3p3i4mKaN28OIJ64d+vWjSZNmhAeHs79+/exsLBg8eLF7N69m8DAQPbt24ejo2OdOQ2Rj0WZTEZ6ejpTp07l8OHDJCUlcfLkSTp27Ag8U6jlfalJkybMnDmTqKgoBg0aRGBgIKNGjSIrK+uF9JIS/5yMjAxGjRrF3r17yc/Pf6GfAzg7O6OhocGZM2f47bffgKrytaKiglWrVrFz505kMhljxoxh/PjxmJqaviCHT506xbx58+jSpQuHDh0iPDychg0bsnbtWn799VcA5s2bR0xMDIsWLeLIkSOcOXOGVatW4e3tTYMGDV54h5omU+oyldf++fn5LF26lGPHjuHi4kJkZCS7du2iSZMmrF27VtzEycvLIzMzU3STunLlCn5+fgwcOFDcwJZbp9YFefquUyt8mB88eICxsfG/KkMumBYtWsTSpUtFJVhfXx99fX00NDS4ffs29+7dw9fXFwUFBVRVVVFTUyMuLo7Lly/TrVu3lwq4ur4z9CraSOL1IrXRu01FRQUPHz4U84AeOnQIQ0ND6tWrh6qqKs7OzhgaGrJixQo+++wz7OzsUFJSIiUlhVatWmFoaIhMJmPjxo00b95cPPGSW7/Udv9l+Wlq5XdMTU2lW7dunDp1ij59+tC/f38cHByws7NjyZIlNG/e/IWgXZUjFzs5OYkRbStTE8aSfK764YcfOHnyJMXFxWhpaaGtrY26ujo3btzg2rVr9OzZE0VFRTGAEcCFCxfIyMjAw8MDJSUl1NTUxAW//LT1Xe9Lr3LdkJqayrx581BSUuL9998nOTmZxMREfH19xbRtlX1bTU1NcXBwoLy8nFu3btGqVSvWrl2Lnp7ev3upWsY/baOKigo0NDTIy8vj/PnzmJiYYG5uXm0wNU1NTbZt24aqqiru7u5VrEru3r1LcHAw7u7uNGvWDFNTU0xMTMSTZ4D4+HiMjIx4/PgxlpaWjBw5EplMRlBQEKGhoWRnZ6OhoYG9vT0mJiZkZWURGhqKgoICkZGRuLq6Mnfu3Jdas7zr1ARZ9zqpHD8EICgoiFatWjFu3DhOnz7N0KFDsbOzo169ejRv3pzt27ejra2Ng4NmCr7XAAAgAElEQVSDGHgzJCSEiIgI9u3bxwcffMCYMWPEuA9y3nV5WheQFOb/JzExkQkTJpCens5nn32GlZUVV65cITo6ml69emFkZERGRgYxMTHUq1dPTFRvaGjI48ePUVJSws3NrUYsFN4GdV2o1gSkNnq3UVBQ4MGDB5w5c4aAgAASExNp2rQpTZs2RVNTE21tbTZv3kxUVBT//e9/mThxIqWlpezduxdtbW1cXFxo2bIlZWVl+Pv7U69evRfKr20cOnSIpk2bitZDz7+jsrIyKioqREZGMmbMGExMTBAEgcaNG3P37l0OHz5M9+7dxYBLcuSKt5GRUZW/5dSEsXTkyBFGjBjBw4cPycnJISoqivj4eLp27YqhoSHFxcVitOtmzZqJUZ/NzMy4cuUKurq6VfIyy6kp/eiftFF1Abf69u3L3r17MTU1Zd68eejp6aGlpcW6deto3Lgxtra2qKioUF5eTklJCUePHsXS0hIzMzPatWtHly5d8PT0BP488Fxd45+0UUVFhZgC1NHRkeDgYHJzc2nRogWamppiG8rr2dbWVhzrd+/eRRAEMjIy+Pnnn/n666+xt7dn6NChYrvLFdtdu3YxevRobty4Qbdu3dDR0cHOzo5bt24REBDAvXv3WLNmDeXl5Rw6dAgbGxuaNm2Kk5MT69evZ9q0aQQEBODj4yM+d00ZO5WpCbLudSLvR5mZmXz++edi2rDGjRsTFRVFjx49sLCwQBAEDA0NSU1NJTo6Gnt7e9q0aUPLli3Jy8tDXV2dpUuX0qZNG6BuBferKdQ5hbm6UwaA69ev8+jRI9avXy8usHbs2MH169epV68ejo6OGBoakpCQQFJSEu7u7qirqyOTyWjZsiWdOnUSU0hIvEhdF6o1AamN3i2eD5QkCALr169n165dTJkyhcmTJ2NsbEy9evXEQFNr1qyhWbNmjBo1CoCIiAhu3LhBRkYGLVq0wMjIiLZt24r31GZ5denSJcaOHYudnR2WlpYA/PLLL+zbt4/09HTs7OxQVlamUaNGREZGUl5ejre3t1gvdnZ2rF+/Hi0tLVq3bl2l7D/zU37Xx9LVq1dZu3YtQ4YMYf78+fTq1YuioiK2bNmChYUFVlZWaGhokJyczJkzZ+jduzdKSkqUlpaioqKCm5sbnp6eNVq5+7vrhsp+yvA/BcfAwIAdO3ZgYWFBz549AVBXVyc7O5u9e/eir6+PpaUlgiAQERFBUFAQtra24m/LT+8rK3ESz/gn40hBQUGMLnz58mWys7O5evUqBgYGtGjRotpTZmdnZ3R1dTl48CBHjx4lPj6ee/fuoaOjw/Dhw7G0tBQ33GJiYhg9ejRnz57F1taWnJwcBgwYgLq6OoqKimzcuJGysjKWL1+OmZkZRUVF4omyjY0NBgYG9OvXD3t7e1RUVMS2r6my+F2Xda+bvLw85s+fz+7duzE2NmbDhg2oqalhZ2dHaGgoWVlZuLm5iabVDg4O7N69m5ycHFq2bImFhQXt2rXDx8enSjqzmtofajN1SjpX9lMuKSmp4p/m7e3NhAkTKC4u5osvvmDkyJG4ubnh7+/PTz/9RGZmJqampnh5eXHz5k1CQ0PFe+UnNZKPgYSExKtCPmnKfeXy8/NFJbB3797Ur18fdXV1cnNzUVBQoLy8nKdPnyKTyUhJSeHWrVukpqbyySefMG3aNDESbOXyazNOTk50796dNWvWiLv/K1eu5N69e8yePZtZs2ZRWFiIqakpo0ePZufOndy9e1dMzdW4cWM++OADvvnmG9LT09/26/xlcnNzxaBEL4uua21tjZ+fH0OGDOHJkydMmzaN1atXY2Fhwbx584Bn6ch8fX158uQJ33//PfC/yMAaGhpA7Zzznn8n+QJWUVGR7Oxsjhw5wo0bN0TltlOnTrRt25bk5GR+//138b7p06fj4uLCt99+y8CBAxk4cCBz5sxh2LBhODk5VfkNSVF+Ncg3NpYvX86HH37IjRs3KC4uJj09nf3794sB+uRtLD8t1tPTY9iwYfzyyy8cOHCAxYsXs3//frp160a7du3Ecr/77juGDRuGr68vx44dw9raGmtra1HxffLkCWfOnKF58+bo6uoCzzburK2tSU1NFf1cGzVqJD6D1PY1h+rih8gtbxISEjAwMEBNTY3i4mIAZs2aRXh4ODExMWJ76+jo0LNnTxISEsTr5MHm3qW0exIvUqdOmOUdceXKlaxcuZKIiAjy8vIwMTFBXV1djPYZHx/PggULGDJkCI8fPyY0NJSKigo8PDwwMzPDzMyMPn361FhTtLdBXd+FrAlIbfTuERISwr59+3BwcEBLS4utW7eSm5tLVlYWwcHB7Nmzh4ULF5KTk4Obmxuqqqr88ssvBAcHs2XLFlxdXQkICMDMzOxtv8pboUWLFnz//ffIZDIUFRVZsWIF/fv3p3Xr1ixatAhzc3PRVPLChQvExMTQo0cP4Jk8b926Nba2ti+cMP8Zb3Ms9erVi19//RUvLy/xVGPv3r3cvXsXbW1tNDQ0UFJSomXLlty/f5+AgACUlZVZvnw57u7u/PzzzygoKNC2bVu0tLR49OgR+vr6ODo6vvBbNXnOe76N5JYFld+p8qny999/z9SpU7l27Rpbt24lLS2NZs2aoa2tjb29PatWrcLS0hJbW1vxHi8vL5ydnWnatCnm5uasWLFC3Liq7RYer4I/G0fyzYzKbZefn8+SJUsICAhg/PjxdO/enebNm7Nr1y50dHRwcHB4qb/wwYMHCQsLw83NDW1tbVq3bo1MJuP27dvo6+tTUVHBZ599hp+fHzKZjO3bt2NlZUXbtm2BZ5YF586d4/z58+Tm5hISEsKvv/7KsmXLGDFihKhEQ80eO5WpC+uG5834T506hUwmQ1lZGS0tLXR0dLhw4QKlpaX4+/uLaWabNGlCXFwcFy9epH379mhqagLg4uLC+++/j7a2dpXfqS19orZSqxXm5yek9PR0xo8fz/Xr1xkxYgTa2tpERERw+fJlOnXqhEwmY9GiRVhbWzNkyBDg2cC4f/8+Fy5cwM/Pj4YNG4rJ6aUJ769TF4RqTUdqo7dLdT5sYWFhnDlzBgMDA6ytrVFVVeXSpUskJibSokULmjdvjpubG9u3b8fa2pru3bvToUMHmjRpwueff0737t2rLChrM9X5f8oXJKtXr8bc3Bx/f38EQcDMzIy7d+9y5MgRvL29MTQ0pFGjRqxevRobGxvRhFZZWVkMjvZ3eFNjqbKLkfz9TUxMWL9+Pa1ataKwsJARI0Zw+vRpzp8/T1hYGM7Ozujr66OgoMBPP/3E/fv3Wbp0KSYmJty/f5+QkBDOnz9Pnz59aNSoEW5ubri4uLz2d3nTPN9G8vERFRXF2bNnadWqFQoKChQWFrJy5Uqio6NZsGAB//nPf/D29iYwMBBVVVVsbGwwNjbm/v37hIWF4enpKSpGysrKmJiY0Lx5c1q3bo2ioqLYTrV9PL4K5G1UnStdaWmpeDJX+fNff/2Vffv2MW7cOAwMDIBnkcmTkpI4e/YsLVq0EAOxPV/upUuX2L17N1ZWVlhbW6OgoMCMGTP4+eefcXJywsHBAQ0NDfGe5cuX07t3bywtLcV82m3btiU1NZWYmBjKyspYsGAB5ubmouVKbWv3urBukLfZ0aNH+eijjzh//jwHDhwgKioKFxcXbGxsKC4u5vjx49jY2GBqair2B7mViZ6enrgBIz9Jro39oTZTa87+K5ueyc0mKvv+AVy7do2nT5+ya9cuevXqRY8ePUhJSeHevXs8ePCAvLw8mjRpwr1797hz5w7Hjx8nNjaWSZMm8csvv4h+cHKkji4hIfF3qSyrSktLxf/LU9M8fPhQ/Ozjjz/GwMCAyMhIUlNTsbS0ZOPGjQQHBzN+/Hg++OAD+vfvj5qamugaYmZmRq9evTA1NaW8vLxOKMsVFRXiqdHt27ermFAPGzYMMzMzysrKKCsrE+eDOXPmkJKSQkhICKWlpbi6uuLl5cWZM2eAd1++V/Z7zcjIICMjg5KSEjp27IijoyM///wzZ86coU+fPhw9epSdO3dSUVHBhg0bRJPt7OxssrOzxejMCQkJfPDBB/Tu3Zvs7GwEQUBVVbVWml4/jyAI7N+/n4CAALZu3SrmKy8rK8PS0pI5c+bg7u7O9evXWbZsGTk5OURERBATEwPA7NmzycrKYu/evX+Y07ymRkN+W1R2pbtx4wZr1qwBnm1GJCUlsWzZMn788Ucxv62NjQ1Pnz4V+3hhYSEA48aN49atW0RERJCdnS2aR8tkMvHaYcOGYWVlRWhoKLdu3QJgwIABZGdni5Hk4Zls+O233ygvL6+SbeDatWvIZDIWLFjAxo0b2bRpExYWFlI+5RrAy1LCyomPj2f58uWMGjWKAwcOEBoayr1795g/fz5PnjyhU6dO2NrasnHjRuBZfygtLcXY2JgpU6bQpk0bcezL5xapP9QsanxrVRZEBQUFrF27lkWLFjFlyhRWr15NWlqa2DkjIyPx9PRETU2NadOm0b17d3x8fFi3bh1NmjRBS0uLrl27kpeXx6hRo5g6dSqenp706NHjBWVZQkJC4u9Q2WetoKCAoKAgduzYIfo9pqSk0LNnTw4cOCAq0lpaWvTq1Yvk5GQiIyPFz0pLS3ny5Am3bt1i/PjxmJiYvJD+CKjVgQjldSQ3m7116xYDBgxgzJgxDBgwgPnz53Pt2jU0NDQICAjgyJEj4oK2vLwcTU1NPvnkE9atW0dCQgLKysoEBgYye/bst/xmfw2ZTMbTp0/56quvGDduHJMnT2bx4sUALFy4kDNnzrB27VratGlDvXr1MDIy4j//+Q+xsbGcO3cOeGYamJ2dzfDhwxkyZAhbtmyhQ4cOfPPNN9jZ2Yl9p7b1oep8uxUUFCgqKkJRURELCwt27NgBgKamJh4eHjg4OLB9+3YxmvqRI0fIysoiPDyctLQ0VFVVGTlyJCEhITx9+vRNv1KtQ95GioqKPH36lM8++4zevXsjCALFxcUsWLCA/v37iyf7X331FTt37kRHR4f33nuPZcuWAf/LY3v37l3q16/P4cOHiYmJQVFRkYKCAj799FMGDBhAYmIiABMnTiQxMZGzZ89SXFyMo6Mjnp6eHDx4kKtXr4rPl56ejpKSEmZmZiQnJzN8+HDGjRtHWloagJhmTfJLfbep3M8KCgrYsmULYWFhYn+QzzOHDh2iRYsWDB8+nPz8fJYvX05mZiYODg6oq6tjbm5O586dSU9PZ8+ePVV+Y8yYMTg7O7/ZF5N45dT4USwXRKtXr8bLy4tLly6hoKBAQUEBGzZsYOLEiZw9exZ4JsB++OEH3N3defz4Mdu2bWPu3LkYGBhw48YNkpKS6Ny5M9u2bWPBggWcOXOGkSNHArUzuImEhMSbQ650bNiwAW9vb86cOUNubi55eXmimXCHDh04ceKEOFnDM59UMzMzjh8/zr1794BnZqNffvklI0aMoEGDBgQFBdGwYcO38l5vmoqKCrZs2cLChQuBZ3PAnTt3+Pzzz7G1tWXNmjWMGzeOK1eu8PXXX5OVlUWfPn2ws7Nj5cqVFBcXizv9Y8eOpVGjRuIplDzP8J+dNrxN5HNReHg4nTt3Ji0tjS+//JIZM2bw8ccfA2BqasrYsWMpKSmpslj38fHB3t6ekJAQfvvtN3x9fVm8eDHm5ubY29sTFhaGu7t7ld+pjcjrJDc3F4CSkhIA3N3dadSoEXp6evz+++8cOnQIAAMDAx4/fsyBAwcYP348s2bNwszMDCsrK3799VeOHTsGPFO2IiMjRVNgib/P88GwfvzxRzw8PCgtLSUsLIwJEyaQlJQk+pIvXbqUH3/8kQYNGrB48WKysrIYNmwYOTk5fPHFF8TFxZGbm8u5c+eYPHkyI0aMwNvbm9WrV+Pm5kZOTg47duzAzs4OAGdnZzw8PKooyJMnTyY3N5ejR4/y5MkTAJKSktDS0mL58uV069aNBg0acPDgQVq2bFnlfWrbZlNto3KMgvbt23P48GHWrFnDsGHDyMjIEAMdZmdno6mpSUhICF27diUpKYndu3czYcIE0brL09MTCwsLNm/eTElJiXgv1G55WmcQajhXrlwR2rdvL/j5+QnHjx8XiouLhZKSEkEQBOHEiRNCnz59BD8/P0EQBOHq1auCl5eXMGPGjCplPHz4UPjiiy+EgwcPivfKKS8vfzMvUsuJjY19248g8SdIbfR6KS0tFb755huhR48eQkREhFBYWCg8ffpUEIT/yZkHDx4IHTt2FJYsWSLk5OSI9546dUpo1aqVMH36dKG8vFx49OiR8PPPPwu3b98Wr6lLsmrWrFlC3759hdOnTwuCIAgHDhwQ2rdvL2RkZIjXHDt2TBgwYICwZMkSQRCe9W87Ozvh8OHDgiAIQllZWZV/XyWveyxlZ2cLgwcPFlauXCmUlpZW+U7+d3l5ueDs7CwsWrRIyM/PF7+/deuW0LFjR2HZsmVCXl5elXvk99V24uLiBD8/P2HcuHFCYWGh+HlCQoLw6aefCiEhIcLYsWOFgIAAsY5CQ0MFT09PIS4uThAEQUhNTRVGjx4tdO7cWVi5cqVQVFQk1t3r6FN1jcjISMHV1VVo3bq1YGNjI6SlpYnfBQYGCh999JEgCIJw/Phxwc/PT/D39xfCwsLEa06ePCn4+fkJXl5eQtu2bYXevXsLjx8/Fq5cuSJ07NhR6Ny5s3DixIlqfzstLU2Uw5mZmYIgCMLGjRsFX19fISIiQhAEQVi3bp1gY2Mj9OnTR7hy5Yp4b10YP5Wp6esGeX/w8/MTTpw4IZSVlQm3b98W3nvvPWH+/PmCIDwbz0FBQULr1q0FLy8vISoqShzjBQUFwpYtW4SkpCRBEJ7JlpSUlLf2PhKvjxp/whwaGkpxcTGzZ8+mY8eOYuQ6eJYqavjw4Tx8+JBNmzZhb29P165dCQkJITg4mLi4OKKjowkICODOnTtYW1tX2RECycdAQkLi1fD48WPOnj3L8OHD8fX1pV69eqirqwOIAUAaNmxI3759iYyM5NKlS+K9GRkZqKurc/HiRS5evIi+vj6DBg3C0tKSioqKOmP2J/z/Lv2QIUPQ0dHhl19+oby8nMePH9OoUaMq8rtt27a0bNmSK1eu8OTJE5ycnPD19WXevHnk5+eLp8zyPLg1iT179nDjxg369esnBj6So6SkRFlZGTKZjClTprBz506uXbsmft+sWTPc3d25du2aeLJa19KaHDp0iAcPHhAVFcWsWbM4deoUAJaWlsTExNC4cWN69+7NgwcPRPNKV1dX8vPz2bRpE1u3bmXSpEnY29uzfft2Jk2ahKqqqlh3kp/yv2PkyJFMmjQJf39/wsPDad++PZ9//jnwzMKksLAQHR0dPvroI/773//Su3dvdu3aRffu3cnJySE/Px9PT0+2bdvGkiVLWLFiBcHBwejp6REaGkpRURGzZ8/G29tb/M179+6xefNmMjMzMTY2pn///kRGRvLrr78C8NFHH6Gnp0dwcDBPnjzB29ubrVu38ssvv9CqVas6JYdrE/L+8NVXX+Ht7Y2ioiKWlpY0atQIZ2dnMW2UPGdy27Zt6dChgzjGz549S1hYGCkpKQC0bt0aU1PTGjenSPw5NT5KtqGhIVeuXCEtLY2OHTuirKwsLqoUFBTQ09Pj1q1bxMXF0atXLzp27EhGRgb79+/n9OnThIWF4efnR2BgYJWQ/xKvlroQSbGmI7XR6+X69ev8+OOPVaK3RkRE8Ouvv3L8+HFyc3OxtLTExcWFsLAw7ty5Q7169VBTU2P//v34+/tjYWFB3759q5T7fCqc2oz8PfX19Xny5Annz5/H0NAQc3NzNmzYQOfOnTE0NERBQQEVFRWSkpKIjY3lww8/RCaT4ezsTMuWLbG2tq623FfF6x5LISEh5OfnM3bsWKD6FEUKCgq0bNmS0NBQUlNTcXV1Ff05vby86N27t/h35XvqAgYGBiQlJWFlZYWWlhYREREoKipib2/P3bt3SU9PZ+DAgcTHxxMfH4+DgwONGzfGwMCA1NRUoqKi6NSpE5999lmVnNR1pf5eN61btyYgIAB9fX0xddemTZuwtLTEysqK9PR01q1bh4WFBUFBQeLar7CwkE2bNvH48WOaN2+OhoYGJiYmmJqaAs/ayNDQkDt37vDbb7/RpUsXAObOncvMmTOxsrLC3d0dFRUVnJ2d2bt3Lzk5Odja2qKtrY2amhoHDhygY8eO2Nra0rhxY7Hcuhr5vKavG6rrD1999RXh4eFkZWVx4MABLC0tcXZ2pqSkhD179nDmzBmSk5PZtm0bP/zwA0OHDqVfv35Vyq2LfaG2U+MVZn19fa5fv05SUhKqqqq0aNECeNZZKyoq0NDQICEhgWvXrtGjRw80NDTw9vamT58+uLi4MH78eDp06ABUn9ZF4tVQ04VqXUBqo9eLiYkJ+/btIyoqitjYWJYsWcL169fFIEzBwcE0aNCAVq1a0bhxY2JiYti3bx/bt2/HyMiITz/9FBUVFTHNSl2VVfJ3NzY25sqVK1y+fJmhQ4cSExNDXFwc7u7uoi9yVFQUxcXFvPfeeygpKaGurk7Tpk1f+zO+7rG0f/9+8ZRLS0urWmVZntbE0tKSwMBAbGxsxFQ58sV9XZ3z9PX1iY+PJyMjA3d3d9zc3Fi2bBllZWVkZmbSoEED3N3dUVBQIC4ujszMTDw8PGjevDm+vr4MGDAALy8voPr8zRL/Dh0dHdTU1EhLS8PExARdXV0ePXrETz/9xKhRo7Czs+PQoUNoa2tjZ2eHnp4eeXl5rF69mlOnTuHj40OTJk1eKFdBQQF9fX0eP35MTEwMFy5cYP78+VRUVLBkyRLef/99VFRUxFNFHR0d9uzZg4GBgbjR1r9//xfGdl1u+5q+bqjcH86fP8/8+fMpLy/niy++wMHBgWPHjhEREYGzszM+Pj5i/vqMjAwaNGjAmjVraNeuHSBtmtV2arzCDM+i2D18+JDExETatWtH/fr1xcipgLhL/PHHH6OkpCSmyjAwMEBVVVVcNEgd/fVR04VqXUBqo9ePh4cHhYWF3Lt3jx49euDv78+YMWMYNmwYv//+O+Hh4YwYMQJTU1N8fHxwcHBg6NChDBs2DCUlJe7fv4+xsXGdllXyvNIaGhqUlZVx7tw5BEHgo48+4rvvviM2Npa8vDwxldLo0aNfCMTzunndY0lZWZkff/wRR0dHLC0tX8i1nZ6eTkBAAB06dMDa2prCwkK8vLzE9FG1Nfr136G0tJSUlBRu3rzJxIkTady4MefOnSM0NJSysjL69euHhYUFCQkJREdH4+TkhIGBgej2Ja0bXj8PHz7E2NgYVVVVDA0N2b9/P4WFhbi6umJtbc358+dZtWoVCQkJrFq1ivT0dObPn4+rq+sflmtiYiK65E2cOJG5c+dWyc0sN7eVp5hSU1PDxcUFRUVFVFRU6uxGU3XUhnVD5f4wYcIE5s+fj4WFBRYWFjg5OREUFET79u2xtLSkcePG+Pj44OPjQ+fOnVFXV5dkQR1B6c8vefdp0KABnTt3ZuvWrezZs4fx48eL3925c4fw8HAGDBhAvXr1qvUxkXxOJCQk3gQWFhZMmTKlWpmjqqpKvXr1xBMuDQ2NFyIWSxPyM+T14OfnR0xMDAcPHqRr165s3ryZHTt2EBERQXl5OevWrauV6TxcXV1p3bo1GzZswNzcHFtb2yp9IzY2Fk1NTTHa97Rp097Wo76z6Orq4ufnxw8//MC2bdv4+OOPad26NdnZ2fj5+Ymb7gMHDqRLly6i9Zocad3wZrGxsWH48OGsXbuWkSNH4uLiQmBgIJcvXyYnJ4c+ffqIJrV/dtJnaGhI165duX//PkVFRcD/0tMpKCgQHx9Pamoq3bt3Z8OGDWKKKDlS29cuKvcHea7t8vJyFBUVKSkpoX79+i/EN1JRUQHqTtwHiVqQVkpO586dad68OadOnRJzbaalpfH9999jZGREjx49AGnBKSEh8XapbnJNTk4mNTWVHj16oKen98I10u71i1RUVKCqqkq3bt2QyWSsW7cOe3t7Fi5cyHfffcfPP/+Ms7OzGIynNqGhocHMmTO5desWM2bMIDo6moSEBBISEpg6dSqBgYH4+flhaGgo3lPb6uBV0LlzZ+zt7Tl+/Djx8fEYGRkRFBTE4MGDxTFoa2uLm5vbW35SCVVVVfz9/bG0tBRzpTdo0ABfX1/69esnKst/9fTXx8eHFi1aEBUVJa4ZHz16xH/+8x8GDBjAzZs3AdDW1hbLlai9PN8fFBUVefToEWvXrqVly5YvtVKS5uW6Q61RmFVVVenevTvKysrs2bOHjRs34u/vT0FBAQsWLKjWn0VCQkLibREXF0d8fDwrVqygf//+NGzYkCFDhrztx6oxyBUaNzc3XF1duXr1KidPngQQT4TKy8trbTCeVq1asXTpUrS0tBg7diyTJk3is88+Izs7m+3bt1cbHE6iKpXXDfv37xc/k5Sjd5PGjRszdOhQDh8+TGpqapU+LTyXv/nPqG7N2LlzZx48eEBERASfffZZlfKkU8Tazcv6Q0lJCd9++63oziJRd6kVJtlyXF1dOXXqFJs2baJhw4asWLFCTBtQ2adZQkJC4m1SVFTE2rVryc7ORlFRkeXLl+Ph4QFIgUP+DvK66tGjBzExMVy7dg1PT0+x/mp7ep8uXbrQpUsXbt++TUFBAfXq1RMjgEtz3l/D1dWVkydPEh4eTkxMDC4uLlK9vaMoKirSuXNnHBwcxMjXcv6JzKxuzVg5CKzUD+oWUn+Q+CNqlcIM0L9/f1xcXERF+e/uOkpISEi8burVq8esWbPIysqiVatWgOSn/E+QB7uytLSktLSU+/fvvxAAqy7QrFkz8f/SnPf36devH25ubri4uLztR5H4E3R0dNDR0XllY1xaM0pURuoPEi+j1inMTZo0EaE1uicAABGySURBVM2vpR0hCQmJdxVTU1PxlESSVf8cBQUFbt26RWlpKTY2NuJndZW6/O7/FHlEXJAsPGoKr6qNpDWjRGWk/iDxMmp1T5A6uoSERE1AklX/jrCwMBwdHRk4cODbfhSJGo6kLNddJDksURmpP0hUptadMEtISEhI1C0mT54sLW4kJCQkJCQkXgvSCkNCQkJCokYjKcsSEhISEhISrwtplSEhISEhISEhISEhISEhUQ2SwiwhISEhISEhISEhISEhUQ2SwiwhISEhISEhISEhISEhUQ2SwiwhISEhISEhISEhISEhUQ2SwiwhISEhISEhISEhISEhUQ2SwiwhISEhISEhISEhISEhUQ2Kc+bMmfO2H+Lv8ujRI65du4aKigr169fnwYMHGBsbv+3HkvgDakobPd+3asNv/9Vyr127RmZm5t/6fXnZubm5JCf/Xzv3H1Nl+f9x/MUBDjAVyxSCrJH9MNcPKzd0pVsZLlDI1Kkrs9K0P5pr1eY0W5mZIv7RZi7/aO1T/8RqbmnTYZtOJ6lT0nC4oLbIkw4UUFOHGOA51/cPd873eHOfc+4Dh3M4nedjYxPPfb+v9329r+s+vIEbT9znraOjQ4cPH9aJEyfU0tKi3Nxc2/H9xzU3NweOCc69vr7+ltes13bgwAHt379fN27cUFtbW+B6Q113uHmxe62rqytQo+B/W2MeOXJEhYWFIefYGjvU+JHmLNKaCRXD7rxQsexqEu5a7OYmOLb/tXBr0Xp88DWcP39eP/74o/7++2/l5eVFzN/62pkzZ1RTU6OsrCz19vb2+35nl2Nzc7N6enps12m480Md09TUpJqaGl27dk1tbW2289nfe4B1/q0xnc6ndY86qanTeG63W1euXAnUyG4dhYvn9H7ndD7D7ZGBvCdEe36kPTkYY1r51+ZA91Esc0p0/FgIntcxY8bELG6yfG1nlQw1wxBjkkx1dbXJyckxI0eONDk5Oaa6utocP3480WkhgmSokd3aSvaxncatrq42WVlZUY3vj52Tk2MkBf4dr3mrrq42mZmZRlLgw+129xm/urrauN3uwDGZmZlmxYoVt+RuPd96bcEf6enpgTh21x1uXuxey8zMNG6324wcOdK43W6TmZnZpw7+84YPHx5yjq2x7fKzzoX/uOB4kdZMqHn3z2nweaFiWWNY62Y3T9a5CR7PP4fh1mJwLm6327hcrj619X+4XK6w+VvjWWMtWLCg32s6+JoyMjJs87Nb505qZ4wxK1as6FN/63z29x7gX1v+/7PGDDeGNXfrHo1UU6fx/J9v2LChz3Gh9p/dtUa630Vz7w21RwbynhDt+ZH25GCMaWVdm/3dR7HMKdHxY8E6rytWrIhZ7GT42s4qGWqGoSepGub29vY+X8Dm5OSYvXv3Jjo1RDDUb6qh1lZ7e3vSju00bn/GtzsnnvPW3t5usrOzI44fLs9QH9nZ2SFjh/vIyckxjY2NYecl2rihYlrn2Ml1hruu7Oxs097eHnEthJt3u5h2sRobG21j+MfpT83CzV9/YwY34P2J19jYGPWajiZHJ2vAekxjY2O/4w8031jVJ9w50cTLysoKu1+dxhvIPTXUcU72fLRrKdz5ofb1YI5pFWptRruPYplTouPHwmDMa7Ch/rWdVTLUDENTUj3D7P8VqGCZmZlqbW1NUEb4rwi1tjweT9KO7TRuf8a3O8fpubHg8XiUnp5u+5rL5QqM7/F45HJFd5tLS0sLGTuczMxM1dXVhZwXl8sVddxQMa1zHK4efunp6UpLSwv5msfjibgWws27XUzr3Puvxy6Gv25OrsUpf+79jWmXv9N4dXV1UY0VbY7B6zzU+dZ1Ek1OA7kHRDtGNLHCnRNNvIyMjLD71Wm8gdxTQx3nZM+HE+09PdS+tq6xWI5pFWptRruPYplTouPHwmDMazJLhpphaMpIdALRKCoqUk9Pzy3/F6tnXJDaQq2toqKipB3badz+jG93jtNzY6GoqEher9f2NZ/PFxi/qKhIPp8vqtjGmJCxw+nt7VVxcXHIefH5fDLGxCSmdY7D1cMv3DV5vd5AvHBjhZt3u5jWBt1/PXYxgusW6VqcCs69PzGtayeaeMXFxVGN5aSG1twirQHrOokmp4HcA/ozhtNY4c6JJt6NGzfC7len8QZyTw11nJM9H0609/RQ+9q6xmI5plWotRntPoplTomOHwuDMa/JLBlqhiEq0T/ijpb/2YPc3NzAswfJ9ishqSgZamS3tpJ9bKdx/c8wRzO+P7b/V/kiPdMXa0PlGWb/rx5bnz20mxdr3Ozs7MDzt7m5uYFnKK118J83bNiwkHNsHXegzzCHWguRnmEOPi9ULKfPMAfPoXVugsfzz2G4tRicS7TPMNvNRfBrsX6G2X9N/X2GOdw+jvQM80DuAf615d8T1pjhxrDmbt2j1n0W6XpDxfN/bn2GOdz+s7vWSPe7aO69ofbIQN4Toj0/ls8w9zfnwXyGebDe1xP5dYNTPMN8q2SoGYaeNGOi/JHHENDR0SGPx6OioiKNGTNGJ06c0KRJkxKdFsJIlhpZ19Z/YWyncfft26eRI0dGNb4/9vDhw9XZ2Rn3eevo6FB9fb0uX76s2267TU888YTt+P7jJAWOCc797Nmzt7xmvbZTp06pra1NTz75pNxud+B6Q113uHmxe01SoEbB/7bG3LNnj8rKykLOsTV2qPEjzVmkNRMqht15oWLZ1STctdjNTXBs/2vh1qL1+OBrGDZsmH799Vfl5+fr2WefjZi/9bULFy6orq5OxcXF6urq6vf9zi5HSbr77rtt12m480Md09TUpLq6Ot1///1yu92289nfe4B1/q0xnc6ndY86qanTeEVFRTpz5kygRnbrKFw8p/c7p/MZbo8M5D0h2vMj7cnBGNPKvzYHuo9imVOi48dC8LxOmDAhZnGT5Ws7q2SoGYaWpGyYrZJ1w6YSajT0UaOhjxolB+o09FGjoY8aDX3UCKkiqf7oFwAAAAAA8ULDDAAAAACADRpmAAAAAABs0DADAAAAAGCDhhkAAAAAABs0zAAAAAAA2KBhBgAAAADABg0zAAAAAAA2aJgBAAAAALBBwwwAAAAAgA0aZgAAAAAAbNAwAwAAAABgg4YZAAAAAAAbNMwAAAAAANigYQYAAAAAwAYNMwAAAAAANmiYAQAAAACwQcMMAAAAAIANGmYAAAAAAGzQMAMAAAAAYIOGGQAAAAAAGzTMAAAAAADYoGEGAAAAAMAGDTMAAAAAADZomAEAAAAAsEHDDAAAAACADRpmAAAAAABs0DADAAAAAGCDhhkAAAAAABs0zAAAAAAA2KBhBgAAAADARpoxxiQ6CQAAAAAAhhp+wgwAAAAAgA0aZgAAAAAAbNAwAwAAAABgg4YZAAAAAAAbNMwAAAAAANigYQYAAAAAwAYNMwAAAAAANmiYAQAAAACwQcMMAAAAAICNpGyYr1+/rnfeeUczZsxQaWmpDhw4YHucz+fTp59+qpkzZ6qiokJvvPGG2tra4pxtanJaI0lqamrSokWLNHPmTM2cOVMHDx6MY6apK5oaSVJ3d7dmzZqluXPnxilDOK3Rvn37NHfuXJWXl2vWrFn63//+F+dMU8/p06e1cOFCPf/881q4cKE8Hk+fY7xer9atW6eSkhLNmDFD27dvj3+iKcxJjb744gvNmjVLFRUVmjt3rn7++ef4J5rCnNTI76+//tLEiRNVVVUVvwThuEY1NTWqqKhQeXm5KioqdOHChfgmCgwmk4S2bt1qPvjgA2OMMadPnzZPPfWU6ezs7HPc3r17zfz5801vb68xxpiNGzeatWvXxjPVlOW0RteuXTPTp0839fX1xhhjent7zaVLl+Kaa6pyWiO/yspK8/7775s5c+bEK8WU57RGJ0+eNOfPnzfGGHP16lVTUlJifvnll7jmmmoWL15sdu7caYwxZufOnWbx4sV9jtmxY4dZunSp8Xq95uLFi2batGnm7Nmz8U41ZTmpUW1trenq6jLGGNPU1GQmTZpkrl+/Htc8U5mTGhljzI0bN8wrr7xi3nvvPbNp06Z4ppjynNSooaHBlJWVmfb2dmPMzfehf//9N655AoMpKX/CvGfPHi1cuFCSVFRUpEceeUS1tbW2x/b09Ki7u1s+n0/Xrl3TnXfeGc9UU5bTGu3evVuTJk3S448/LknKyMjQ7bffHtdcU1U0++j48ePyeDyaPXt2PFNMeU5rNHHiROXn50uSRowYofvuu08tLS1xzTWVXLx4UY2NjSovL5cklZeXq7GxUZcuXbrluJqaGs2fP18ul0ujRo1SSUmJfvrpp0SknHKc1mjatGnKycmRJI0fP17GGF2+fDnu+aYipzWSpC+//FLPPPOMioqK4pxlanNao2+++UZLly7VmDFjJN18H8rKyop7vsBgScqGubW1VXfddVfg84KCAp0/f77PcdOnT1dxcbGmTp2qp59+WqdPn9bSpUvjmWrKclqjP//8UxkZGVq+fLlmz56tNWvW6MqVK/FMNWU5rVFXV5c2btyodevWxTM9yHmNgjU3N+vkyZOaMmXKYKeXss6dO6f8/Hylp6dLktLT05WXl6dz5871Oa6wsDDwuZP6ITac1ijYzp07dc899/CN9ThxWqPff/9dhw4d0uuvv56ALFOb0xo1Nzfr7NmzWrRokebMmaNt27bJGJOIlIFBkZHoBOzMmTNHra2ttq8dOXLEcZzffvtNzc3Nqq2t1bBhw7RhwwZt2rRJH330UaxSTVmxqpHP59PRo0f13XffafTo0aqsrNSmTZtUWVkZq1RTVqxqtHnzZr388svKz88P+3wZoherGvm1t7frrbfe0tq1awM/cQYQWV1dnbZs2cLz/0NMb2+vPvzwQ1VWVgaaNgw9Xq9Xf/zxh77++mv19PRo2bJlKiws1Isvvpjo1ICYGJIN844dO8K+XlhYqJaWFo0aNUrSze+ATZ482TbOlClTNGLECEnSCy+8oDVr1sQ+4RQUqxoVFBRo8uTJysvLkyRVVFRQoxiJVY1OnDih2tpabdu2Td3d3bpy5YoqKiq0a9euQck7lcSqRtLNX51bsmSJli1bprKyspjniv9XUFCgtrY2eb1epaeny+v1qr29XQUFBX2Oa21t1WOPPSap70+cMXic1kiS6uvrtXLlSm3btk3jxo1LQLapyUmNOjo6dObMGb355puSpKtXr8oYo87OTq1fvz5RqacMp/uosLBQpaWlcrvdcrvdeu6559TQ0EDDjP+MpPyV7NLSUn3//feSJI/Ho1OnTmnatGl9jhs7dqyOHj2q3t5eSdLBgwf1wAMPxDXXVOW0RmVlZWpoaFBnZ6ckqba2VuPHj49rrqnKaY127dql/fv3a//+/frss8/04IMP0izHidMa/fPPP1qyZIkWLVqk+fPnxzvNlHPHHXdowoQJ2r17t6Sbf4thwoQJgW9s+JWWlmr79u3y+Xy6dOmS9u3bp+effz4RKaccpzVqaGjQu+++q88//1wPP/xwIlJNWU5qVFhYqGPHjgXeg1577TUtWLCAZjlOnO6j8vJyHTp0SMYY9fb26ujRo3rooYcSkTIwKNJMEj5k0NXVpdWrV6upqUkul0srV65USUmJJGnLli3Ky8vTSy+9pO7ubn388cc6efKkMjIyVFBQoPXr1/OrinHgtEbSzefGvvrqK6WlpWns2LFav369Ro8encj0U0I0NfI7duyYqqqq9MMPPyQi5ZTjtEZVVVX69ttvde+99wbOffXVVzVv3rxEpf6f19zcrNWrV+vq1avKzc1VVVWVxo0bp+XLl+vtt9/Wo48+Kq/Xq08++USHDx+WJC1fvjzwR9ww+JzUaN68eWppabnl64LNmzfzjds4cVKjYFu3blVXV5dWrVqVoIxTj5Ma+Xw+VVVVqba2Vi6XS1OnTtWqVavkciXlz+WAPpKyYQYAAAAAYLDxrR8AAAAAAGzQMAMAAAAAYIOGGQAAAAAAGzTMAAAAAADYoGEGAAAAAMAGDTMAAAAAADZomAEAAAAAsEHDDAAAAACAjf8DEl/vwFfACucAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 864x72 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "selected_authors = np.array(\n",
    "    [\n",
    "        \"Dean Heller (R)\",\n",
    "        \"Bernard Sanders (I)\",\n",
    "        \"Elizabeth Warren (D)\",\n",
    "        \"Charles Schumer (D)\",\n",
    "        \"Susan Collins (R)\",\n",
    "        \"Marco Rubio (R)\",\n",
    "        \"John Mccain (R)\",\n",
    "        \"Ted Cruz (R)\",\n",
    "    ]\n",
    ")\n",
    "\n",
    "sns.set(style=\"whitegrid\")\n",
    "fig = plt.figure(figsize=(12, 1))\n",
    "ax = plt.axes([0, 0, 1, 1], frameon=False)\n",
    "for index in range(authors.shape[0]):\n",
    "    ax.scatter(authors[\"ideal_point\"][index], 0, c=\"black\", s=20)\n",
    "    if authors[\"name\"][index] in selected_authors:\n",
    "        ax.annotate(\n",
    "            author_map[index],\n",
    "            xy=(authors[\"ideal_point\"][index], 0.0),\n",
    "            xytext=(authors[\"ideal_point\"][index], 0),\n",
    "            rotation=30,\n",
    "            size=14,\n",
    "        )\n",
    "ax.set_yticks([])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "qnVQbe7l2gQG"
   },
   "source": [
    "## Automatic guide generation\n",
    "\n",
    "Above, for pedagogical reasons, we defined the guide (i.e. Variational Family) manually, making sure it is identical to that descibed in the original paper. \n",
    "\n",
    "However, manually defining the guide is often unnecessary as NumPyro contains a [module](https://num.pyro.ai/en/stable/autoguide.html#id8) that enables automatic guide generation based on the model provided.\n",
    "\n",
    "In our case it turns out that `AutoNormal` creates a guide that is identical to the one we manually defined above. Specifically it first transform variables to belong to unrestricted space of real numbers. For example, by applying the log transformation to those variables that are restricted to be non-negative (document locations $\\theta_d$ and objective topic locations $\\beta_j$). Then it uses independent Normal distribution as a guide for each transformed distribution.\n",
    "\n",
    "In the cell below we verify that the guide generated by `AutoNormal` is in fact identical to the guide manually defined above as a part of the TBIP class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "id": "DvUOdaVKu-1k"
   },
   "outputs": [],
   "source": [
    "from numpyro.infer.autoguide import AutoNormal\n",
    "\n",
    "\n",
    "def create_svi_object(guide):\n",
    "    SVI(\n",
    "        model=tbip.model,\n",
    "        guide=guide,\n",
    "        optim=adam(exponential_decay(learning_rate, num_steps, decay_rate)),\n",
    "        loss=TraceMeanField_ELBO(),\n",
    "    )\n",
    "\n",
    "    Y_batch, I_batch, D_batch = tbip.get_batch(\n",
    "        random.PRNGKey(1), counts, author_indices\n",
    "    )\n",
    "\n",
    "    svi_state = svi_batch.init(\n",
    "        random.PRNGKey(0), Y_batch=Y_batch, d_batch=D_batch, i_batch=I_batch\n",
    "    )\n",
    "\n",
    "    return svi_state\n",
    "\n",
    "\n",
    "# This state uses the guide defined manually above\n",
    "svi_state_manualguide = create_svi_object(guide=tbip.guide)\n",
    "\n",
    "# Now let's create this object but using AutoNormal guide. We just need to ensure that\n",
    "# parameters are initialized as above.\n",
    "autoguide = AutoNormal(\n",
    "    model=tbip.model,\n",
    "    init_loc_fn={\"beta\": initial_objective_topic_loc, \"theta\": initial_document_loc},\n",
    ")\n",
    "svi_state_autoguide = create_svi_object(guide=autoguide)\n",
    "\n",
    "\n",
    "# Assert that the keys in the optimizer states are identical\n",
    "assert svi_state_manualguide[0][1][0].keys() == svi_state_autoguide[0][1][0].keys()\n",
    "\n",
    "# Assert that the values in the optimizer states are identical\n",
    "for key in svi_state_manualguide[0][1][0].keys():\n",
    "    assert jnp.all(\n",
    "        svi_state_manualguide[0][1][0][key] == svi_state_autoguide[0][1][0][key]\n",
    "    )"
   ]
  }
 ],
 "metadata": {
  "accelerator": "GPU",
  "colab": {
   "collapsed_sections": [],
   "provenance": []
  },
  "gpuClass": "premium",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.13"
  },
  "vscode": {
   "interpreter": {
    "hash": "72bfc691e125b99acbc6d80c2190d4963aecde9ca99ebcb34dfd45ffa5b3f8fc"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
